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FOREWORD 


The  developments  documented  in  this  report  were  carried  out 
at  the  Aeroelastic  and  Structures  Resea. ch  Laboratory,  Department 
of  Aeronautics  and  Astronautics,  Massachusetts  Institute  of 
Technology,  Cambridge,  Massachusetts  under  Contract  *Io.  FM620-70- 
C-0020  from  the  Office  of  Scientific  Research,  U.3.  Air  Force. 

The  software  system  originated  in  connection  with  Dr.  Orringer's 
need  to  use  finite-element  analysis  for  problems  in  the  compressive 
behavior  of  fiber  composites  that  were  being  investigated  under  the 
Contract.  Dr.  Jacob  Pome^antz  of  the  Aeromechanics  Division,  AFOSR 
served  as  technical  monitor. 

The  authors  wish  to  express  their  appreciation  to  Prof.  Pin 
Tong  of  the  Department  of  Aeronautics  and  Astronautics  for  the 
numerous  helpful  suggestions  he  made  regarding  the  design  of  the 
FEABL  software  system.  The  authors  are  also  indebted  to  all  of  the 
students  in  the  Departments  1971-72  finite  element  analysis  course 
who,  under  Prof.  Tong's  guidance,  acted  as  guinea  pigs  while  the 
software  was  still  being  tested.  For  his  review  of  this  report  and 
for  attendant  constructive  suggestions,  the  authors  are  indebted  to 
Professor  Emmett  A.  Witmer  of  the  MIT-ASRL. 

The  computations  and  testing  of  the  program  were  carried  out 
at  the  MIT  Information  Processing  Center. 
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FEABL,  a  basic  software  system  developed  for  finite  element 
analysis  at  the  MIT  Aeroelastic  and  Structures  Research 
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not  limited  to)  the  specialized  type  of  continuum  analysis 
problem  encountered  in  the  materials  laboratory.  FEABL  has 
also  been  used  successfully  as  an  educational  tool  in  the 
finite  element  analysis  course  given  by  the  MIT  Department  of 
Aeronautics  and  Astronautics.  The  software  is  modular,  and  is 
written  in  machine-independent  F0RTRAN  IV.  A  complete  program 
listing  is  contained  in  Appendix  C.  Element  subroutines  which 
can  be  used  either  independently  or  in  conjunction  with  FEABL 
will  be  presented  in  future  publications. 
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SECTION  1 


INTRODUCTION 


1.1  Purpose  and  Scope 

Finite  element  analysis  programs  developed  in  the  past 
fall  into  two  general  categories:  the  "one-shot"  program 
optimized  for  solution  of  a  particular  problem,  or  the  "systems" 
program  designed  for  the  non-optimum  solution  of  any  problem. 

The  structural  engineer  who  has  a  new  analysis  to  perform  must 
spend  a  great  amount  of  time  either  writing  his  own  "one-shot"  pro¬ 
gram  or  learning  how  to  use  one  of  the  general  "systems".  FEABI 
attempts  to  relieve  the  analyst  of  these  burdens  by  providing 
a  basic  library  which  handles  the  numerical  operations  common 
to  all  finite  element  analysis,  and  which  can  be  understood 
and  used  after  a  few  hours  of  study.  FEABL's  current  capa¬ 
bilities  are  limited  to  static  analysis.  Other  options  will 
be  added  in  the  future,  as  the  need  arises,  and  according 
to  their  usefulness.  This  guide  is  written  for  the  analyst 
who  has  some  acquaintance  with  the  theory  of  finite  element 
analysis  and  with  IBM  FORTRAN  IV. 

THE  FEABL  software  system  has  been  designed  primarily  for 
"materials  laboratory"  stress  analysis.  The  analysis  problems 
encountered  in  this  application  area  are  characterized  by  geometry 
which,  although  not  simple,  is  mathematically  describable,  and 
by  stress  solution  accuracy  requirements  which  can  usually  be  met 
with  models  of  1,000  or  fewer  degrees  of  freedom.  Stress  concen¬ 
tration  analysis  for  two-material  systems  and  for  composite-to- 
metal  step-lap  joints  are  two  examples  of  the  materials  laboratory 
type  of  problem.  Another  feature  of  many  of  these  problems  is  a 
requirement  for  nonstandard  solution  methods.  For  instance, 
the  boundary  conditions  may  include  an  undetermined  prescribed 
displacement,  the  value  of  which  depends  upon  satisfaction  of  an 
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auxiliary  condition  on  a  boundary  stress  integral;  a  specific 
example  is  given  in  Ref.  1.  The  main  features  of  FEABL,  all¬ 
in-core  computation  and  control  of  the  solution  process  by 
the  user,  are  well  adapted  to  the  needs  of  such  problems. 

Experience  in  the  MIT  Department  of  Aeronautics  and 
Astronautics  has  also  shown  that  FEABL  is  a  valuable  educational 
tool.  A  budget  of  $50  to  $70  per  man  on  an  IBM  370/155  machine 
allows  each  student  to  solve  one  or  two  problems,  each  having 
200  to  300  degrees  of  freedom.  In  the  approach  taken  at  MIT, 
the  student  is  required  to  program  his  structure  geometry,  element 
interconnections,  and  subroutines  for  generating  stiffnesses  and 
stresses.  He  must  also  achieve  a  general  understanding  of  the 
principles  upon  which  the  FEABL  subroutines  are  based.  Thus,  each 
student  is  able  to  gain  some  practical  experience,  as  well  as 
theoretical  knowledge  in  one  semester  of  a  finite  element  analysis 
course.  This  approach  is  felt  to  be  superior  to  the  traditional 
method  of  having  the  student  learn  only  how  to  code  the  input  data 
for  a  particular,  complete  "systems”  program. 

The  flexibility  ana  economy  of  the  FEABL  software  system 
has  been  amply  demonstrated  by  various  applications  in  which  real 
problems  were  solved  during  the  test/exercise  development  phase. 
The  following  are  examples  abstracted  from  these  applications: 


Problem 

Total 

CPU 

Total 

Description 

DOF 

Time  (Min) 

Cost  ($) 

Displacement  analysis 

216 

0.42 

5.37 

of  two  large  rigid  frame 
structures  (6  load 

426 

2.49 

13.17 

cases  each) 

Continuum  stress 
concentration  analysis 

338 

1.50 

4.50 

with  boundary  condition 
determined  by  auxiliary 
stress  integral  (3 
complete  solutions 
required  for  each  case) 

442 

2.70 

12.00 

Stress  analysis  of 

1,266 

3.0 

14.00 

a  large  frame  in  an 
oil  tanker,  using 
rectangle  and  triangle 
continuum  elements 
and  flange  elements 
(6  loading  cases) 

(Approx. ) 

!  i  l 

Pressure  distribution  160  I  No  data  available 

analysis  for  viscous  flow 

around  a  corpuscle  I 

in  a  blood  vessel 


Stress  analysis, 
including  solution 
for  stress  intensity 
factor  of  a  sharp 
crack,  using  hybrid 
■element  at  the  crack 

160 

0.13 

2.54 

The  above  examples  were  run  on  the  IBM  370/155  at  the  MIT  Infor¬ 
mation  Processing  Center.  The  costs  presented  are  for  production 
runs,  exclusive  of  debugging. 
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1.2  General  System  Concepts 

Development  of  the  FEABL  software  has  been  based  on  three 
general  concepts.  First,  F0RTRAN  IV  was  selected  as  the  program 
language,  and  all  F0RTRAN  syntax  rules  were  followed  rigorously. 
This  makes  FEABL  as  machine-independent  as  possible,  to  allow 
for  use  at  installations  having  other  than  IBM  equipment.  The 
only  restrictions  are  that  the  F0RTRAN  compiler  software  must 
permit  logical  IF  statements,  variable  dimensioning,  six-character 
variable  names,  and  labelled  C0MM0N  declarations.  Second,  the 
software  is  modular  in  character.  The  finite  element  analysis 
process  has  been  decomposed  into  a  series  of  distinct  steps, 
and  each  step  which  is  independent  of  structure  geometry  has  been 
programmed  as  a  separate  subroutine.  FEABL  is  thus  analogous, 
in  a  sense,  to  the  IBM  Scientific  Subroutine  Package. 

The  third  concept  deals  with  data  storage  techniques.  In 
order  to  minimize  wasted  storage  space,  a  vector  approach  has 
been  adopted.  All  of  the  input  data  for  a  problem,  certain 
"housekeeping"  data,  and  internally  generated  problem  data  such 
as  the  master  stiffness  matrix  are  stored  in  a  single  vector, 
referred  to  below  as  the  /DATA/  vector.  In  addition,  control 
parameters  which  locate  entries  in  the  /DATA/  vector  or  manage 
input-output  operations  are  organized  in  four  C0MM0N  groups. 

1.3  General  System  Organization 

The  analyst  must  provide  one  or  more  programs  which  inter¬ 
face  with  the  FEABL  software: 

1.  A  MAIN  program  which  inputs  the  required  problem  data 
and  controls  the  execution  of  all  subprograms. 

2.  One  or  more  subroutines  which  generate  an  element 
stiffness  matrix  (and,  if  required,  equivalent  nodal 
forces  for  thermal  stress,  gravity,  etc.)  for  each 
type  of  element.* 


Subroutines  from  ASRL  EGL  (Element  Generator  Library)  may  be 
used. 


3.  One  or  more  subroutines  which  generate  element  stresses 
from  the  element’s  nodal  displacements.* 

Data  may  be  input  by  reading  cards,  disk  or  tape  files,  by  auto¬ 
matic  generation  techniques,  or  by  a  combination  of  the  two 
methods.  The  details  of  a  particular  problem  will  determine 
which  approach  is  more  efficient.  For  very  simple  types  of  prob¬ 
lems,  in  which  the  structure  is  divided  into  only  a  few  different 
element  types  and  shapes,  stiffness  matrix  and  stress  generation 
may  be  done  within  the  MAIN  program.  However,  if  separate  sub¬ 
routines  are  required  for  these  tasks,  pre-existing  ones  may  be 
interfaced  with  FEABL,  since  the  FEABL  software  package  has  been 
written  in  a  manner  such  that  pre-existing  F0RTRAN  subroutines 
may  be  adapted  to  run  with  FEABL  with  little  or  no  re-programming. 

The  interface  between  user-written  programs  and  FEABL 
consists  of  two  parts:  data  location  and  process  control. 

Figure  1  illustrates  the  general  features  of  the  data  location 
interface.  The  labelled  C0MM0N  areas  mentioned  in  the  previous 
section  must  appear  in  the  user's  MAIN  program.  This  is  accom¬ 
plished  simply  by  including  identical  sets  of  DIMENSI0N,  C0MM0N 
and  EQUIVALENCE  declarations  in  MAIN  and  in  the  FEABL  software. 
There  is  no  requirement  for  the  element  generator  subroutines 
to  communicate  with  FEABL  in  this  manner;  these  subroutines  need 
only  be  interfaced  with  user’s  MAIN  by  an  identical  DIMENSI0N 
declaration  for  the  element  stiffness  matrix  and  "force/displace¬ 
ment"  vector.  The  latter  is  a  dual-purpose  vector  in  which 
element  equivalent  nodal  forces  are  generated,  and  which  serves 
later  as  a  storage  area  for  the  solution  displacement  vector  of 
the  element.  The  interface  between  the  generator  routines  and 
FEABL  is  achieved  automatically  by  variable  DIMENSI0N  declarations 
contained  in  the  FEAEL  software.  The  data  location  interface  is 
discussed  in  detail  in  Section  2. 

Subroutines  from  ASRL  EGL  (Element  Generator  Library)  may  be 

used. 


5 


User 

/  1/i.l  llall  W  A  1^1  \ 

J  for  clement  stiffness  \ 

Element 

MAIN 

* \  matrix  and  force- 

\olisp  lace  went  vector/ 

generator 

Subroutines 

Labelled  C0MM0N  Areas: 

•  /DATA/  Vector 

•  Con  fro)  Parameters 


FIGURE  1 


The  process  control  Interface  is  simply  a  recognition 
that  input  of  problem  data  and  operations  on  the  data  must  be 
performed  in  a  certain  sequence.  This  requires  a  definite  series 
of  programming  steps  in  user's  MAIN  program.  Section  3  discusses 
the  process  control  interface  in  detail. 
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SECTION  2 


THE  FEABL  DATA  STORAGE  SYSTEM 


2,1  Accessory  Data 

The  accessory  data  is  divided  into  four  control  parameter 
C0MM0N  groups.  The  C0MM0N  declaration  statements  listed  below 
for  the  four  control  parameter  groups  must  appear  in  all  user- 
written  MAIN  programs. 

2.1.1  Input-Output  Control  Parameters 

C0MM0N  /I0/  KR,  KW,  KP,  KT1,  KT2,  KT3 

The  six  parameters  in  the  /I0/  group  are  hardware  device 
codes,  defined  as  follows: 

KR  -  card  reader 

KW  -  printer 

KP  -  card  punch 

KT1,  KT2 ,  KT3  are  extra  device  codes  which  the  user  may  define 
optionally  as  he  pleases.  For  example,  KT1  might  refer  to  a 
system  direct  access  disk,  KT2  to  an  external  tape  unit,  etc. 

Only  the  printer  device  code  KW  is  used  directly  by  the  FEABL 
software  presented  in  this  guide.  The  input-output  control 
parameters  are  used  in  F0RTRAN  READ  and  WRITE  statements,  e.g.: 

READ  ( KR , 1 )  NET,  NDT 

1  F0RMAT  (216) 

WRITE  ( KW , 2 )  NET,  NDT 

2  F0RMAT  (16HOT0TAL  ELEMENTS® ,  16,  10X,  1OHT0TAL  D0F=,  16) 

The  correct  values  of  KW  and  any  other  device  codes  employed  by 
the  user  must  be  established  at  the  beginning  of  the  MAIN  program. 

2.1.2  Problem  Size  Parameters 

C0MM0N  /SIZE/  NET,  NDT 

These  parameters  define  the  problem  size  in  terms  of  the 
total  number  of  elements  (NET)  and  the  total  number  of  degrees 
of  freedom  (NDT)  in  the  whole  structure.  NDT  Includes  both 
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constrained  and  unconstrained  degrees.  In  standard  solutions 
these  parameters  are  established  at  the  beginning  of  the  MAIN 
program  and  remain  fixed. 

2.1.3  Address  Indexi  Parameters 

C0MM0N  /BEGIN/  IC0N,  IK0UNT,  ILNZ,  IMASTR,  IQ,  IK 

C0MM0N  /END/  LC0N,  LK0UNT ,  LLNZ,  LMASTR,  LQ,  LK 

These  parameters  control  the  begin  and  end  locations  of 
data  sub-blocks  in  the  /DATA/  vector.  The  user  establishes  cor¬ 
rect  values  for  these  parameters  by  calling  PEABL  subroutine 
SETUP  after  he  has  established  the  /SIZE/  block.  Use  of  the 
address  index  parameters  to  locate  entries  in  the  /DATA/  vector 
is  discussed  in  Subsection  2.2. 

2.2  The  /DATA/  Vector 

DIMENSI0N  REAL(xxxx),  INTGR(xxxx) 

C0MM0N  /DATA/  REAL 

EQUIVALENCE  (REAL  (1),  INTGR(l)) 

The  above  declarations  serve  to  define  the  /DATA/  vector 
as  a  one-dimensional  array  occupying  a  C0MM0N  area  labelled 
/DATA/  and  having  two  reference  names:  REAL  for  floating  point 
entries  and  INTGR  for  integer  entries.  These  declarations  must 
appear  in  the  user-written  MAIN  program.  The  user  chooses  the 
dimension  integer  "xxxx"  to  suit  his  particular  needs.  The 
DIMENSI0N  declaration  with  user’s  value  for  the  length  of  the 
/DATA/  vector  must  be  inserted  in  each  FEABL  subprogram,  immedi¬ 
ately  following  the  subroutine  name  declaration,  e.g.: 

SUBR0UTINE  ASEMBL  (LNUM,  NDE,  ELK,  ELQ) 

DIMENSI0N  REAL(IOOO),  INTGR(IOOO) 

Some  compilers  restrict  the  length  of  any  vector  to 
32,768  words  (8,000^Ej,)  .  If  this  restriction  is  encountered, 
and  the  user  desires  a  longer  /DATA/  vector,  say  50,000  words, 
the  following  form  of  the  declaration  statements  can  be  employed 

DIMENSI0N  REAL( 25000),  INTGRC 25000 ) ,  DUMMY( 25000) 

C0MM0N  /DATA/  REAL,  DUMMY 

EQUIVALENCE  (REAL(l),  INTGR ( 1) ) 
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The  C0MM0N  /DATA/  declaration  must  then  be  changed  to  the  above 
form  in  each  FEABL  subprogram. 

The  dual  nature  of  the  /DATA/  vector  must  be  recognized 
clearly.  For  example,  the  100th  entry  in  the  vector  may  be 
treated  as  either  a  floating  point  or  an  integer  quantity  in 
arithmetic  or  logical  instructions  by  operating,  respectively, 
with  REAL(IOO)  or  INTGR(IOO).  This  property  will  lead  to  com¬ 
pilation  errors  if  the  user  attempts  to  save  storage  space  by 
declaring  INTGR  to  be  an  array  of  2-BYTE  words  (half-words), 
while  REAL  is  left  as  an  array  of  normal  4-BYTE  words  (full  words), 
or  if  the  user  attempts  to  run  in  double  precision  mode  for 
floating  point  arithmetic. 

The  /DATA/  vector  is  organized  into  six  blocks,  the  limits 
of  which  are  determined  by  the  values  stored  in  the  address 
index  parameters  (Subsection  2.1.3).  In  normal  usage,  the  first 
four  blocks  contain  integer  data  and  the  last  two  contain  float¬ 
ing  point  data.  Storage  conventions  are  detailed  in  the  follow¬ 
ing  sections,  in  the  order  in  which  the  blocks  appear  in  the 
/DATA/  vector. 

2.2.1  Constraint  Vector 

This  is  an  integer  block  which  contains  the  global  number 
of  each  degree  of  freedom  at  which  a  displacement  it;  to  be 
prescribed.  Entries  in  the  constraint  vector  are  referred  to’ 
by: 

INTGR(I)  where  IC0N<I<LC0N 

The  following  example  illustrates  user  action  involving  the  con¬ 
straint  vector 

• 

• 

INTGR ( IC0N) =2 

INTGR( IC0N+1) =25 

INTGR(LC0N)=6 
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By  means  of  the  above  instructions  the  user  has  specified  that 
global  displacements  2,  25,  and  6  will  be  prescribed.  The  pre¬ 
scribed  displacements  need  not  be  listed  in  any  particular  order. 
Any  excess  space  in  the  constraint  vector  is  filled  with  zeros  by 
FKABL. 

2.2.2  Address  Count  Vector 

This  is  an  integer  block  which  contains  information  re¬ 
lating  to  the  absolute  address  in  the  /DATA/  vector  of  each 
diagonal  entry  of  the  master  (global)  stiffness  matrix.  Entries 
in  the  address  count  vector  are  referred  to  by: 

INTGR(I)  where  IK0UNT<I<LK0UNT 

To  obtain  the  relevant  data  for  the  Nth  row,  the  correct  subscript 
is: 

J=IK0UNT+N-1 

INTGR(J)  contains  address  information  for  row  N 

Normally,  the  user  has  no  direct  communication  with  the 
address  count  vector;  its  contents  are  generated  internally  by 
FEABL.  The  absolute  address  of  each  diagonal  entry  is  output 
for  use  in  debugging.  After  the  output,  the  data  in  the  address 
count  vector  are  modified  as  follows: 

Contents  =  (Absolute  Address  of  K^)  -  i 

Subsequently,  any  stored*  entry  can  be  obtained  by 
referring  to: 

REAL(KADR) 

where  the  address  is  given  by: 

KAPP.  =  INTGR(  TK0UNT+I-1)  -hJ 

2.2.3  Leading  Non-Zero  Entry  Vector 

The  LNZ  vector  is  an  integer  block  containing  the  number  j 
of  the  column  in  which  appears  the  leading  non-zero  entry 


Many  entries  of  the  master  stiffness  matrix  are  not  stored  in 
the  /DATA/  vector.  See  Subsection  2.2.6  for  details. 
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of  the  master  stiffness  matrix  for  each  row  i.  Entries  in  the 
LNZ  vector  are  referred  to  by: 

INTGR(I)  where  ILNZ<I<LLNZ 
INTGR( ILNZ+N-1)  refers  to  row  N. 

Normally,  the  user  has  no  direct  communication  with  the 
LNZ  vector;  its  contents  are  generated  internally  by  FEABL,  and 
are  output  for  use  in  debugging  and  evaluating  displacement 
numbering  strategies. 

2.2.4  Master  Assembly  List 

This  is  an  integer  block  containing  all  of  the  information 
which  relates  the  user’s  local  degree  of  freedom  numbers  to  the 
user’s  global  numbering  system.  The  master  assembly  list  is,  in 
effect,  a  Boolean  logical  transformation  between  the  element 
displacement  vectors,  which  together  form  a  linearly  dependent 
set,  and  the  global  displacement  vector,  which  is  an  independent 
set.  Entries  in  the  master  assembly  list  are  referred  to  by: 

INTGR(I)  where  IMASTR<I<LMASTR 

The  master  assembly  list  need  not  be  filled,  but  if  it  is  not, 
a  zero  must  be  stored  immediately  following  the  last  active 
position. 

The  master  assembly  list  is  subdivided  into  two  sections, 
as  indicated  in  Figure  2. 


Location  ir>  /DATA/  Vector a 
IMASTR  XMASTR+NET 


Zero  must  be 
Stored  by  ustr 
in  first  excess  LMASTR 
location  I 


Point  trs 


(One  storage 
location  per  element) 


Global 
DOF  Numbers 


Excess 


Storage 
(if  any  exists) 


FIGURE  2 


11 


The  pointer  section  consists  of  one  location  per  element,  from 
INTGR(IMASTR)  to  TNTGR(IMASrR+NET-l) ,  with  the  convention  that 
these  locations  correspond  to  the  user’s  elements  in  ascending; 
order  1,  2,...,  NET.  Each  pointer  contains  the  absolute  address 
in  the  /DATA/  vector  of  the  location  where  the  assembly  list 
for  an  element  starts,  as  indicated  in  Figure  3. 
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1 


Element  i 


Global  Members 


Element  Z 


Global  Numbers' 


Element  3 


Global  (lumbers 


Pointers  *for: 
Element  3 
Element  7L 
Element  i 


FIGURE  3 


The  remainder  of  the  master  list  contains  the  user's 
element-by-element  sequence  of  global  displacement  numbers. 

As  an  illustrative  example,  consider  the  set  of  elements  shown 
in  Figure  4,  consisting  of  two  ^I-node  rectangles  and  two  3-node 
triangles,  each  with  two  degrees  of  freedom  per  node.  The  total 
length  required  for  the  master  assembly  list  is  given  by: 


Pointers:  1  per  element - M 

Rectangles:  8  DOF/element - 16 

Triangles:  6  DOF/element - 12 


TOTAL - - - 32 


Suppose  that  IMASTR=101.  Then  the  assembly  information,  correctly 
stored,  would  appear  as  shown  in  Figure  5. 
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The  following  crude  set  of  F0RTRAN  instructions  might  be  used  to 
store  the  assembly  information  shown  above: 

NEXT=IMASTR+4 

INTGR(IMASTR)=NEXT  (establishes  pointer  for  1st  element) 
INTGR(NEXT)=1 

INTGR(NEXT+1) =2 

• 

• 

INTGR ( NEXT+7 ) =4 
NEXT=NEXT+8 

INTGR ( IMASTR+1) =NEXT  (establishes  pointer  for  2nd  element) 
INTGR (NEXT) =5 

INTGR(NEXT+1)=6 

« 

• 

INTGR( NEXT+5 ) =8 
NEXT=NEXT+6 

INTGR ( IMASTR+2) =NEXT  (establishes  pointer  for  3rd  element) 

• 

• 

The  use  of  pointers  in  the  master  list  allows  assemblies 
involving  as  many  different  types  of  elements  with  different 
total  numbers  of  degrees  of  freedom  as  desired.  Also  the  use 
of  DOF  numbers  rather  than  node  numbers  permits  assembly  of 
elements  having  different  numbers  of  displacements  at  various 
nodes  without  requiring  any  special  programming  or  conventions. 
The  user  may  choose  any  element  and  displacement  global  numbering 
schemes  and  any  local  displacement  numbering  conventions  he 
desires.  The  only  restrictions  are  that: 

1.  There  must  be  no  gaps  in  the  master  assembly  list. 

2.  The  array  must  be  filled  or,  if  excess  storage  exists, 
a  zero  must  be  stored  after  the  last  displacement 
numb  er . 

3.  All  element  numbers  and  degree  of  freedom  numbers 
must  be  positive. 

4.  The  lowest  element  number  and  the  lowest  global  degree 
of  freedom  number  must  each  be  unity. 
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2.2 . 5  Force/Displacement  Vector 

This  is  a  floating  point  block  which  contains  ail  of  the 
force  and  displacement  information  required  for  analysis  of  a 
structure.  Entries  in  the  force/displacement  vector  are  referred 
to  by: 

REAL(I)  where  IQ<I<LQ 

The  Nth  entry  in  the  vector  is  the  force  or  displacement  associ¬ 
ated  with  the  Nth  global  degree  of  freedom. 

Various  types  of  information  are  overlayed  in  the  force/ 
displacement  vector.  First,  if  the  structure  being  analyzed  is 
loaded  by  continuum  body  forces  (e.g.,  gravity)  or  is  in  a  thermal 
environment,  the  resulting  element  equivalent  nodal  forces  must 
be  assembled  along  with  the  element  stiffness  matrices.  FEABL 
subroutine  ASEMBL  uses  the  force/displacement  vector  as  a  storage 
area  for  this  purpose.  Second,  the  user  must  introduce  his  glo¬ 
bal  concentrated  nodal  forces  and  prescribed  displacements  into 
the  force/displacement  vector.  Finally,  the  FEABL  solution  sub¬ 
programs  store  the  displacement  solution  in  the  force/displac- 
ment  vector  by  overwriting  the  prescribed  quantities. 

The  following  simple  algorithm  enables  the  user  to  com¬ 
municate  with  the  Nth  global  degree  of  freedom: 

NN-IQ+N-1 

REAL(NN)= . . . 

or 

. . ,=f (REAL(NN) ) 

For  example,  suppose  the  structure  in  Figure  4  is  to  be  given 
the  boundary  conditions  shown  in  Figure  6.  Global  degrees  of 
freedom  4  and  8  have  concentrated  forces  A  and  B  applied,  respec¬ 
tively,  while  degrees  12  and  15  have  displacements  C  and  D 
prescribed,  respectively.  Displacements,  1,  2,  5»  6,  13,  and  14 
are  prescribed  to  be  zero.  The  following  set  of  FORTRAN 
instructions  specify  these  boundary  conditions: 
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(a)  Establishment  of  the  constraint  vector 

INTGR (IC0N)*1 
INTGR(IC0N+1)=2 
INTGR(IC0N+2}=5 
INTGR ( IC0N+3) =6 
INTGR(IC0N+J|)=12 
INTGR ( IC0N+5 ) =13 
INTGR  ( IC0N+6 )  =l*i 
INTGR (IC0N+7) =15 

(b)  Input  of  the  prescribed  displacements: 

REAL ( IQ) =0. 

R£AL( IQ+1) =0 . 

REAL( IQ+4) =0 . 

• 

• 

REAL(  IQ+13)  =0 . 

REAL( IQ+11) =C 
REAL (IQ+1 4) =D 

(c)  Input  cf  the  prescribed  forces: 

REAL( IQ+3) =REAL( IQ+3)+A 
•REAL ( IQ+7 ) =REAL( IQ+7 ) +B 

The  instructions  (b)  and  (c)  have  been  written  assuming  that 
element  equivalent  nodal  forces  have  been  assembled  into  the 
force/displacement  vector.  Thus,  all  degrees  at  which  displace¬ 
ments  are  prescribed  must  be  set  to  their  correct  values,  while 
nonzero  concentrated  global  forces  must  be  added  to  the  pre¬ 
existing  assembled  element  forces. 


C,D 


are 

are 


concert  trated  forces 
prescribed  displacements 


FIGURE  6 
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2.2.6  Master  Stiffness  Matrix 


This  floating  point  block  contains  the  essential  entries 
of  the  master  stiffness  matrix,  stored  one  row  after  the  other. 
Entries  in  the  master  stiffness  matrix  are  referred  to  by: 

REAL(I)  where  IK<I<LK 

Normally,  the  user  is  not  required  to  communicate  directly  with 
the  master  stiffness  matrix,  but  an  understanding  of  its  detailed 
organization  will  be  helpful  in  debugging. 

Since  a  stiffness  matrix  is  always  symmetric,  only  its 
lower  triangle  need  be  kept  in  storage.  Thus,  the  general  or¬ 
ganization  of  the  master  stiffness  array  in  the  /DATA/  vector 
may  be  represented  by  the  diagram  shown  in  Figure  7. 
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FIGURE  7 

Signifies,  ^t  additional  savings  in  storage  may  be  realized  for 
problems  with  many  degrees  of  freedom  by  taking  advantage  of 
the  fact  that  a  master  stiffness  matrix  is  normally  banded  and 
sparsely  populated.  The  boundary  line  of  the  shaded  area  in 
Figure  8  represents  the  leading  non-zero  entry  locations  in  a 
hypothetical  stiffness  matrix,  FEABL  subprograms  which  operate 
on  the  master  stiffness  matrix  incorporate  logic  instructions 
which  cause  the  operation  to  be  skipped  if  the  entry  lies 
in  the  unshaded  area  of  Figure  8.  Thus,  the  leading  zero  entries 
for  each  row  are  not  stored,  and  this  is  where  the  address  count 
vector  and  the  LNZ  vector  come  into  play.  The  actual  organiza¬ 
tion  of  the  master  stiffness  array  consists  of  a  sequence  of 


17 


FIGURE  8 

variable-length  sections  for  the  rows  of  the  stiffness  matrix. 
Figure  9  illustrates  a  sample  sequence. 


Column 


The  algorithms  for  operating  on  the  master  stiffness  matrix, 
based  on  the  contents  of  the  address  count  and  LNZ  vectors  are 
quite  simple: 

(For  an  operation  on  K . . :  row  I,  column  J,  J<I) 

•  ij  — 

• 

INDEX=ILNZ+I-1 

IF  (J  .LT,  INTGR( INDEX))  G0  T0  5 

KADR=IK0UNT+I-1 

KADR=INTGR ( KADR) +J 

(To  define  address  of  in  master  stiffness  array) 

(Operate  using  REAL(KADR)  for  ) 

5  (Skip  operation  if  J<LNZ  column  number  for  the  row) 

2,3  Modification  of  the  FEABL  Data  Storage  System 

The  data  storage  system  presented  above  has  been  designed 
for  '’production''  computing.  It  has  been  assumed  implicitly  that 
the  user  will  be  conducting  a  great  number  of  studies  involving 
similar  stress  analysis  problems  and  employing  nearly  the  same 
number  of  degrees  of  freedom.  Thus,  the  /DATA/  vector  need  be 
DIMENSI0Ned  only  once,  after  which  production  object  decks  of 
the  FEABL  software  may  be  made. 

However,  the  FEABL  software  may  also  be  placed  in  on-line 
storage  in  a  form  which  will  handle  problems  of  widely  varying 
size,  with  only  minor  modifications.  These  modifications  are 
as  follows: 

1.  Delete  the  C0MM0N  /DATA/  REAL  declaration  from  all 
programs  and  subprograms . 

2.  Delete  the  EQUIVALENCE  (REAL(l),  INTGR(l))  declaration 
from  all  subprograms.  (Retain  this  declaration  in 
MAIN.) 

3.  Use  the  standard  declaration: 

DIMENSI0N  REAL( 2 ) ,  INTGR(2) 
in  all  FEABL  subroutines.  ( DIMENS I0N  the  /DATA/ 
vector  properly  in  MAIN.) 

^ .  Add  the  array  names  REAL  and  INTGR  as  arguments  of 


all  FEABL  subroutines,  e.g.: 

SUBR0UTINE  0RK( LENGTH,  REAL,  INTGR) 

SUBR0UTINE  FACTPD(REAL,  INTGR) 

When  the  on-line  version  of  FEABL  is  used,  only  the  DIMENSI0N 
declaration  for  the  /DATA/  vector  in  the  user’s  MAIN  program 
need  be  changed  to  perform  analyses  requiring  different  amounts 
of  data  storage. 
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SECTION  3 


FEABL  PROCESS  CONTROL 


3.1  General 

For  present  purposes  the  finite  element  analysis  of  a 
structure  will  be  divided  into  eight  programming  stages: 

1.  Establishment  of  input/output  device  codes,  problem 
size,  and  address  index  parameters. 

2.  Input  of  the  master  assembly  list  and  organization 
of  the  master  stiffness  matrix  into  corresponding 
segments . 

3.  Generation  and  assembly  of  element  stiffnesses  (and 
nodal  equivalent  forces,  if  any)  in  global  coordinates. 

ij.  Application  of  rotation  transformations  to  the  master 
stiffness  matrix  (and  assembled  nodal  equivalent  forces) 
at  any  nodes  at  which  the  boundary  conditions  are 
to  be  given  in  special  coordinate  systems. 

5.  Input  of  the  prescribed  quantities,  i.e.,  global  num¬ 
bers  at  which  displacements  are  to  be  prescribed, 
values  of  prescribed  displacements,  and  accumulation 
of  values  of  any  nonzero  global  concentrated  forces. 
Application  of  constraints  to  the  master  stiffness 
matrix  and  force/displacement  vector. 

6.  Solution  for  the  master  displacement  vector. 

7.  Application  of  inverse  rotation  transformations  to 
the  master  displacement  vector,  at  any  nodes  where 

a  rotation  was  applied  in  stage  to  produce  a  master 
displacement  vector  entirely  in  the  global  coordinate 
system. 

8.  Extraction  of  element  displacement  vectors  from  the 
global  vector;  calculation  of  element  stresses  from 
the  element  displacement  vector. 
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Specific  FEABL  subprograms  are  associated  with  each  of  the  above 
stages,  according  to  the  following  table: 


Stage 

FEABL  Subroutines 

5 

BC0N 

6 

FACTPD/FACTSD,  SIMULQ 

hh 

R0TATE 

XTRACT 

Stage 

FEABL  Subroutines 

1 

SETUP 

2 

0RK 

3 

ASEMBL 

H 

R0TATE 

Figure  10  illustrates  the  standard  FEABL  process  sequence  in  terms 
of  the  eight  stages  described  above. 

Analysts  who  are  just  beginning  to  work  with  FEABL  are 
advised  to  observe  strictly  the  standard  process  sequence  out¬ 
lined  above.  This  sequence  will  be  described  in  full  detail 
in  Section  4.  More  experienced  analysts  may  find  it  convenient 
to  depart  occasionally  from  the  standard  process  in  order  to 
reduce  program  length. 

3.2  Description  of  the  FEABL  Software 

The  FEABL  package  consists  of  the  eight  subroutines  listed 
in  the  table  in  Subsection  3.1.  Each  subroutine  is  associated 
with  a  particular  stage  in  the  standard  process  sequence  (except 
R0TATE ,  which  is  associated  with  both  stages  4  and  7).  All  FEABL 
subprograms  are  ready  to  use  for  in-core  finite  element  analysis, 
once  the  DIMENSI0N  declaration  for  the  /DATA/  vector  has  been 
inserted  (see  Subsections  2.2  and  2.3).  The  following  subsections 
demonstrate  how  each  of  the  FEABL  subroutines  is  called,  define 
the  subroutine  arguments,  and  outline  briefly  what  the  subroutine 
does.  A  summary  table  of  the  C0MM0N  area  information  required 
by  each  FEABL  subroutine  is  given  in  Appendix  A.  Listings  appear 
in  Appendix  C. 
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FIGURE  10 


3.2.1  Housekeeping  Setup  Subroutine  for  /DATA/  Vector 
(SETUP) 


CALL  SETUP( LENGTH,  NC0N,  MASTRL) 

LENGTH  -  A  scalar  integer  numerically  equal  to  the  dimen¬ 
sion  which  the  user  has  assigned  to  the  /DATA/ 
vector. 

NCON  -  A  scalar  integer  greater  than  or  equal  to  the 

total  number  of  degrees  of  freedom,  in  the  assem¬ 
bled  structure,  at  which  displacements  are  to 
be  prescribed. 

MASTRL  -  A  scalar  integer  greater  than  or  equal  to  the 
total  number  of  words  required  for  the  master 
assembly  list. 

Based  on  NC0N,  MASTRL  and  the  total  number  of  degrees 
of  freedom  in  the  entire  assembled  structure  (NDT,  in  the  /SIZE/ 
group) ,  subroutine  SETUP  organizes  the  /DATA/  vector  by  calcu¬ 
lating  the  address  index  parameters  in  the  /BEGIN/  and  /END/ 
groups,  except  for  the  index  LK  which  defines  the  end  of  the 
block  reserved  for  the  master  stiffness  matrix.  SETUP  uses  the 
argument  LENGTH  to  test  whether  the  user’s  /DATA/  vector  has 
at  least  enough  storage  available  to  accommodate  the  first  four 
data  blocks  (constraint  vector,  address  count  vector,  LNZ  vector 
and  master  assembly  list).  If  the  /DATA/  vector  is  too  short, 
SETUP  estimates  the  total  length  required  for  all  six  data  blocks, 
based  on  a  reasonable  population  density  for  the  lower  triangle 
of  the  master  stiffness  matrix,  prints  the  estimate  and  aborts 
the  run.  If  the  first  four  blocks  can  be  accommodated,  the  con¬ 
straint  vector  is  filled  with  zeros  and  control  is  returned  to 
MAIN. 
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3.2.2  Subroutine  for  Detailed  Organization  of  the  Master 
Stiffness  Matrix  Block  (0RK) 

CALL  ORK( LENGTH) 

LENGTH  -  A  scalar  integer  numerically  equal  to  the  dimen¬ 
sion  which  the  user  has  assigned  to  the  /DATA/ 
vector. 

Subroutine  0RK  produces,  in  essence,  a  map  for  the  master 
stiffness  matrix  like  that  of  Figure  9,  using  the  information 
contained  in  the  master  assembly  list  to  calculate  the  correct 
values  of  the  entries  in  the  LNZ  vector.  This  is  accomplished 
by  first  setting  each  LNZ  column  number  equal  to  its  row  number 
(diagonal  matrix),  and  then  examining  the  assembly  information 
element  by  element  to  re-set  the  LNZ  column  numbers,  according 
to  the  following  algorithm: 

1.  The  smallest  global  number  N  for  the  element  is  found. 

2.  The  element’s  global  numbers  are  then  treated  as 
row  numbers.  If  the  LNZ  column  number  corresponding 
to  a  row  (global  number)  is  greater  than  N,  its  value 
is  re-set  to  N. 

Once  the  LNZ  vector  has  been  established,  subroutine  0RK  uses 
an  accumulation  process  to  calculate  the  absolute  address  of 
the  diagonal  entry  of  each  row.  By  convention,  the  diagonal 
is  the  only  entry  stored  for  the  first  row  and  is  therefore  located 
by  the  address  index  parameter  IK: 

INTGR( IK0UNT) =IK 

Subsequent  entries  are  located  by  the  algorithm: 

INTGR( IK0UNT+M-1 ) =INTGR( IK0UNT+M-2 ) +M+1-INTGR ( ILNZ+M-1 ) 

(Diag  Addr  for  Mth  row)=(Diag  Aadr  for  row  M-l)+(Total  no. 

of  nonzero  entries  in  Mth  row) 

Conveniently,  the  address  of  the  diagonal  for  the  last 
row  in  the  master  stiffness  matrix  is  also  the  correct  value 
of  LK.  Subroutine  0RK  now  tests  the  /DATA/  vector  by  means  of 
the  argument  LENGTH.  If  the  /DATA/  vector  is  too  short  for  the 
problem  data,  an  exact  calculation  of  its  required  length  is 
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output  and  the  run  is  aborted.  If  sufficient  storage  is  avail¬ 
able,  the  stiffness  matrix  map  is  output  and  the  address  count 
vector  is  modified  by: 

INTGR ( IK0UNT+M-1 ) =INTGR ( IK0UNT+M-1 ) -M ;  1 <M<NDT 
for  all  rows,  M.  This  saves  repeated  subtraction  of  the  row 
number  in  later  subprograms.  The  correct  algorithm  for  locating 
Kmn  in  the  /DATA/  vector  is  now: 

KADR=INTGR( IK0UNT+M-D+N;  M>N 

K,„  is  assigned  to  REAL(KADR) 

Since  the  next  stage  of  the  analysis  will  involve  the  accumulation 
of  data  in  the  master  stiffness  matrix  (and  perhaps  in  the  force/ 
displacement  vector),  0RK’s  last  action  before  returning  control 
to  MAIN  is  to  fill  these  two  data  blocks  with  floating  point 
zeros . 

3.2.3  Element  Assembly  Subroutine  ( ASEMBL) 

CALL  ASEMBL (LNUM,  NDE,  ELK,  ELQ) 

LNUM  -  A  positive  scalar  integer  =  user’s  global  element 
number 

NDE  -  A  scalar  integer  =  total  number  of  degrees  of  free¬ 
dom  possessed  by  the  element  which  is  about  to 
be  assembled 

ELK  -  A  floating  point,  two-dimensional  square  array 

which  contains  the  stiffness  matrix  of  the  element 
about  to  be  assembled. 

ELQ  -  A  floating  point  vector  which  contains  the  equiva¬ 
lent  nodal  forces  for  the  element  about  to  be  assem¬ 
bled,  or  which  contains  floating  point  zeros  if 
there  are  no  equivalent  nodal  forces. 

Since  ELK  and  ELQ  are  variably  DIMENSI0Ned  in  this  sub¬ 
program,  elements  having  different  numbers  of  degrees  of  freedom 
can  be  handled  automatically.  However,  these  arguments  must 
be  DIMENSI0Ned  explicitly  in  the  user’s  MAIN  program.  For  example, 
suppose  a  structure  is  to  be  analyzed  in  plane  stress  with  a 
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combination  of  3-node  triangle  elements  (6  degrees  of  freedom) 
and  4-node  rectangle  elements  (8  degrees  of  freedom).  Then  the 
declaration: 

DIMENS I 0N  TRIK(6,6),  TRIQ(6),  RECK(8,8),  RECQ(8) 
might  appear  in  MAIN.  A  triangle  element  would  be  assembled 
by  the  instruction: 

CALL  ASEMBL (LNUM,  6,  TRIK,  TRIQ) 
while  a  rectangle  element  would  be  assembled  by: 

CALL  -ASEMBL  (LNUM,  8,  RECK,  RECQ) 

Subroutine  ASEMBL  can  handle  elements  having  as  many  as  one  hun¬ 
dred  degrees  of  freedom  (not  a  serious  restriction).  If  the 
assembly  of  larger  elements  is  attempted,  ASEMBL  will  abort  the 
run  and  tell  the  user  to  change  the  DIMENSI0W  of  one  of  its 
internal  parameters. 

ASEMBL  examines  the  section  of  the  master  assembly  list 
belonging  to  element  number  LNUM  and  records  the  values  of  the 
global  displacement  numbers  in  an  internal  vector  called  MNUM. 
Then,  each  entry  ^  of  the  lower  triangle  of  the  element  stiff¬ 
ness  matrix  (i>J),  is  accumulated  into  its  proper  place  in  the 
master  stiffness  matrix,  according  to  the  algorithm: 

I-MNUM(i) 

J=MNUM( j ) 

*'or  ^ 

Kij"*KJI  for  I<J 

(See  Subsection  3.2.1  for  address  algorithm  for  K^j.)  ASEMBL 
does  not  use  the  upper  triangle  (i<j)  of  the  element  stiffness 
matrix;  the  user  may  omit  calculating  these  entries  to  save  time. 
The  vector  ELQ  of  element  equivalent  nodal  forces  is  accumulated 
in  the  same  manner  into  the  proper  locations  in  the  master  force 
vector  (fifth  block  of  the  /DATA/  vector). 


27 


3.2.4  Rotation  Transformation  Subroutine  (R0TATE) 


CALL  R0TATE(N0DE,  IR0W,  JR0W,  KR0W,  Z ANGLE,  Y ANGLE,  XANGLE) 
N0DE  -  A  scalar  integer  "convenience"  number  for  the  user. 
May  be  positive  or  negative  but  not  zero.  (See 
explanation  below.) 

IR0W  ►  Three  scalar  integers  equal  to  the  global  numbers 
JR0W  of  the  degrees  -of  freedom  at  the  node  at  which 
KR0WJ  the  rotation  is  to  be  performed. 

ZANGLE\.  Floating  point  values  of  the  three  Euler  rotation 
YANGLE  angles  in  degree  measure. 

XANGLE , 

The  argument  N0DE  does  not  enter  directly  into  the  trans¬ 
formation  calculations,  but  is  printed  out  in  the  heading  of 
the  information  supplied  by  R0TATE.  Normally,  the  user  assigns 
either  his  node  numbers  or  the  series  1,  2,  3,...  to  a  set  of 
rotations.  A  positive  value  of  N0DE  causes  the  full  subprogram 
to  execute,  and  is  used  for  stage  4  in  the  process  sequence. 

In  stage  7,  only  the  solution  displacement  vector  requires  trans¬ 
formation;  a  negative  value  of  N0DE  will  cause  execution  of  the 
now  unwanted  transformation  of  the  stiffness  matrix  to  be  skipped. 

The  global  numbers  of  the  degrees  of  freedom  at  the  node 
being  rotated  may  be  contiguous  (e.g.,  13,  14,  15)  or  separated 
(e.g.,  5»  10,  15)  depending  upon  the  global  numbering  scheme 
the  user  has  adopted.  Either  form  is  acceptable  to  this  subpro¬ 
gram.  If  only  two  degrees  of  freedom  are  to  participate  in  the 
rotation  (as,  for  example,  in  plane  stress  problems),  an  integer 
zero  must  be  specified  for  the  third  global  number.  If  six  de¬ 
grees  of  freedom  are  to  participate  (e.g.,  shell  elements),  two 
separate  rotations  are  required,  i,e.,  one  CALL  for  the  trans¬ 
lational  and  one  for  the  rotational  degrees  of  freedom.  Errors 
in  the  global  number  arguments  (e.g.,  repeating  a  number  or  too 
many  zeros)  will  cause  an  abort. 


The  three  Euler  angle  arguments  must  be  specified  according 
to  the  conventions  illustrated  in  Figure  11.  Let  XYZ  be  the 
user’s  global  cartesian  axis  system,  with  respect  to  which  the 
element  stiffness  and  equivalent  nodal  forces  have  been  generated 
and  assembled.  Let  XYZ  by  a  local  coordinate  system  with  respect 
to  which  the  displacement  constraints  are  to  be  specified. 


FIGURE  11 


The  following  conventions  have  been  adopted  in  subroutine  R0TATE: 

1.  Axes  X,  Y  are  first  rotated  through  angle  ZANGLE  about 
axis  Z  to  the  intermediate  positions  x,  y. 

2.  Axes  x,  Z  are  then  rotated  through  angle  YANGLE  about 
axis  y,  x  to  the  final  position "X  and  Z  to  an  inter¬ 
mediate  position  z. 

3.  Axes  y,  z  are  then  rotated  through  angle  XANGLE  about 
axis  X  to  their  final  positions  Y,  Z. 

4.  Positive  angles  obey  the  "right  hand  rule"  of  vector 
analysis . 

Subroutine  ROTATE  forms  the  matrix  of  direction  cosines: 

cos(X,  X)  cos(X,  Y)  cos(X,  Z) 

cos(Y,  X)  cos(Y,  Y)  cos(Y,  Z) 

cos(Z,  X)  cos(Z,  Y">  cos(Z,  Z) 


from  the  values  of  the  Euler  angles  supplied  in  the  argument 
list.  Then,  if  argument  N0DE  is  positive,  R0TATE  applies  the 
following  transformations  to  portions  of  the  master  stiffness 
matrix  K  where  i,  j,  k  are  the  global  numbers  specified 
by  the  user: 
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Finally,  the  corresponding  entries  of  the  force/displacement 
vector  are  transformed  according  to: 


A  few  examples  will  illustrate  the  proper  use  of  subroutine 
R0TATE.  First,  suppose  that  global  degrees  of  freedom  1  and  2 
are  to  be  rotated  by  *J5  degrees  at  node  1  in  a  plane  elasticity 
problem.  The  correct  instruction  is: 

CALL  R0TATE(1,  1,  2,  0,  H5.,  0.,  0.) 

To  return  the  force/displacement  vector  to  the  user’s  global 
axis  system  after  the  displacement  solution  has  been  obtained: 

CALL  R0TATE(-1,  1,  2,  0,  0.,  0.) 

Rote  that  only  the  first  angular  argument  ZANGLE  is  used  in  two- 
dimensional  problems.  In  a  three-dimensional  problem,  to  rotate 
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degrees  4,  5,  6  at  node  2  through  angles  of  15,  30  and  45  degrees 

a/ 

about  the  Z,  y  and  X  axes: 

CALL  R0TATE(2,  4,  5,  6,  15.,  30.,  45.) 
and  after  the  displacement  solution  has  been  obtained: 

CALL  R0TATE(-2S  4,  5,  6,  -45.,  -30.,  -15.) 

Note  that  not  only  the  signs,  but  also  the  order  of  application 
of  the  angles  is  reversed.  However,  some  care  is  required  in 
three  dimensions.  Suppose  degrees  7,  8,  9  at  node  3  were  rotated 
only  about  axes  Z  and  y: 

CALL  R0TATE(3,  7,  8,  9,  20.,  40.,  0.) 

Then  to  reverse  the  rotation  after  obtaining  the  displacement 
solution: 

CALL  R0TATE(-3,  7,  8,  9,  -40.,  -20.,  0.) 

Note  that  the  unused  XANGLE  does  not  participate  in  the  reversal. 
The  instruction: 

CALL  R0TATE(-3,  7,  8,  9,  0.,  -40.,  -20.) 
would  be  incorrect. 

3.2.5  Boundary  Constraint  Subroutine  (BC0N) 

CALL  BC0N 

^  * 

Let  U}  U  represent  respectively  the  global  displacement 

vector  and  the  global  force  vector.  The  result  of  the  first 
four  program  stages  has  been  to  supply  a  right-hand  side 
and  the  stiffness  coefficients  for  the  force-displacement 
relations : 

Km  c  d 

However,  the  value  of  u  is  known  at  some  degrees  of  freedom, 
with  the  corresponding  Q  unknown,  while  the  value  of  Q  is  known 


The  term  "global”  should  be  taken  in  a  more  general^sense  here, 
It  refers  to  the  final  set  of  coordinate  systems  y?Zz  in  which 


the  user  will  prescribe  his  boundary  conditions. 
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at  the  other  degrees. 

Conceptually,  the  set  of  all  degrees  of  freedom  in  the 
structure  may  be  divided  into  two  subsets:  Those  at  which  forces 
are  prescribed  (P)  and  those  at  which  displacements  are 
prescribed  (D).  When  the  user  specifies  the  values  of  the  pre¬ 
scribed  quantities,  Q  becomes  a  "force/displacement"  vector 
in  fact: 


V  A  means  a. 
prescribed  quantity) 


and: 


Kp  =  unknown  displacements 
Q ^  =  unknown  reaction  forces 


The  force-displacement  relations  may  be  partitioned  in  a  similar 
manner: 


K 


FD 


W 


Subroutine  BC0N  transforms  the  force-displacement  relations  from 
the  above  form  to: 


With  the  above  "right  hand  side"  in  the  force/displacement  vector, 
standard  equation-solving  techniques  may  be  used  to  produce  the 
solution  displacement  vector.  Before  making  the  transformation, 

BC0N  arranges  the  entries  of  the  constraint  vector  (Subsection  2.2.1) 
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In  ascending  order  and  checks  for  the  presence  of  gJLobal  numbers 
(positive  integers).  The  run  is  aborted  if  no  global  numbers 
are  found. 

3.2,6  Stiffness  Matrix  Factoring  Subroutine  (FACTPD/FACTSD) 

CALL  FACTPD 
or 

CALL  FACTSD 

The  force-displacement  relations  are  solved  by  Choleski's 
direct  method,  which  consists  of  two  programming  steps: 

1.  The  master  stiffness  matrix  is  factored  into  a  triple 
product. 

2.  The  displacements  are  solved  for  sequentially,  in 
three  sub-steps. 

Subroutine  FACTPD/FACTSD  accomplishes  the  first  programming  step. 
The  master  stiffness  matrix  K  is  factored  into  the  form: 

K-LDlT 

where  L.  is  a  lower  triangular  matrix: 


l 

0 

0 

0  4  4  4 

L21 

l 

0 

0 ... 

L  ~ 

L31 

L32 

1 

0.  • . 

L4l 

• 

L42 

• 

L43 

• 

1  •  •  • 

• 

• 

V 

• 

• 

• 

and  D  =  f D1  D2...  )  is  a  diagonal  matrix.  The  factor¬ 
ing  algorithms  are: 
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d.  11  * 

^  ^*3Yhj»)  t  v  r 


im-1 

2 
^  -  vT  (m) 


Im=  Kw„-  S 


U'  i;Z^  ...^  NbT 


where : 

J(n)  3  Leading  Nonzero  Entry  Column  Number  for  Row  n 
J(m)  =  Leading  Nonzero  Entry  Column  Number  for  Row  m 
J(m,n)  =  max(J(n),  J(m)) 

The  master  stiffness  matrix  is  destroyed  and  replaced  by  the 

entries  of  L  and  D  as  the  factoring  process  is  executed. 

The  entries  L  (m>n)  and  D  are  stored  respectively  at  K  (m>n) 
mn  m  mn 

and  K  .  The  unit  diagonal  entries  of  are  not  stored, 
mm 

FACTPD/FACTSD  tests  the  entries  of  D  for  nonsingularity 
and  positive-definiteness  as  they  are  created.  K  is  positive- 
definite  if  all  D  >0.  Rows  for  which  D  <0  are  reported.  If  any 
Dm  is  found  to  equal  zero  exactly,  K  is  singular;  the  row  in  which 
the  singularity  was  discovered  is  reported  and  the  run  is  aborted. 

The  names  FACTFD  and  FACTSD  refer  to  different  entry  points 
in  this  subprogram.  If  K  is  supposed  to  be  positive-definite 
(as  in  the  case  of  a  structure  analyzed  by  compatible  displace¬ 
ment  elements),  FACTPD  should  be  called.  If  errors  were  made 
in  the  assembly  or  constraint  stages  cf  the  program,  they  will 
appear  now  as  Dm<0  in  one  or  more  rows,  and  FACTPD  will  abort 
the  run.  If  K  is  not  necessarily  positive-definite  (as  in 
the  case  of  a  structure  analyzed  by  hybrid  stress-displacement 
elements),  FACTSD  should  be  called.  FACTSD  continues  execution  even 

if  there  exist  D  <0.  If  the  FEABL  software  is  to  be  converted 

m 

to  the  on-line  storage  version  (Section  2.3),  be  sure  to  make 
the  array  names  REAL,  INTGR  arguments  of  both  entry  points. 
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FACTPD/FACTSD  also  makes  a  rough  estimate  of  the  condition 
ing  of  K  by  calculating  the  so-called  rounding  error  parameter 
(see  Ref.  2,  pg.  8l ) : 

E  =  n>in  (  ,/lKwml) 

E  is  a  measure  of  how  many  significant  figures  of  information 
have  been  lost  in  the  diagonal  entries,  as  a  result  of  the  factor 
ing  algorithm.  The  run  is  aborted  if  E<10  .  The  user  must 

keep  in  mind  the  fact  that  the  rounding  error  parameter  is  an 
imperfect  conditioning  measure.  If  FACTPD/FACTSD  reports  that 
no  significant  figures  have  been  lost,  it  does  not  necessarily 
follow  that  there  is  no  error  in  the  displacement  solution. 
Advanced  techniques  for  realistic  solution  error  estimates  are 
discussed  in  Section  5. 

3.2.7  Subroutine  for  Solution  of  Simultaneous  Equations 
(SIMULQ) 

CALL  SIMULQ (ENERGY) 

ENERGY  -  A  scalar  floating  point  variable,  the  value  of 
which  is  undefined  when  SIMULQ  is  called. 

With  the  stiffness  matrix  in  factored  form  K* LDL , 
let  P =  LT Vi  and  Rc  BP  .  Then: 

UR  e  Q 

can  be  solved  sequentially  for  R1,  R2,...,  RNDT  via  the  algorithm 

RrQi 

Rm  =  Qvn  (  Lywj  Rj  )  }  Wl  e  Z,3,  .,NDT 


Then: 


>vi  e  1,2,  .. .  j  NDT 
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Finally,  the  displacement  solution  is  obtained  by  solving 
sequentially  for  uNDT,  uNDT-1,...,  u1  according  to: 


^not  “ 


U  -  P 

1  v. 


'to 


to 


NOT 

-  r  (l 


V* 


to*  blDT-lj  hiDT-2 ,  . i 


As  SIMULQ  carries  out  these  three  sub-steps,  the  prescribed 
vector  Q  is  first  replaced  by  R  ,  then  R  is  replaced  by  P 
and  finally  P  is  replaced  by  u  .  Q  and  u  are  printed 
out  by  SIMULQ. 

Since  the  approximate  value  of  the  strain  energy  in  the 
structure  is  often  useful  to  the  analyst,  this  quantity  is  cal¬ 
culated  by  SIMULQ  during  execution  of  the  solution  steps.  If 
the  factored  form  of  K  is  introduced  into  the  strain  energy 
expression,  there  results: 

Strain  Energy=  ~  uTKu  |PtDP 

Thus,  the  straightforward  algorithm: 

MDT 

Strain  Energy^  £  Z!  T>:  (Pj) 

ri  3  3 

can  be  used.  The  value  of  the  strain  energy  thus  calculated 
is  printed  out.  At  the  end  of  execution  of  subroutine  SIMULQ, 
the  force/displacement  vector  (fifth  block  in  the  /DATA/  vector) 
contains  the  master  displacement  solution  vector,  and  the  strain 
energy  value  has  been  assigned  to  the  argument  ENERGY. 

3.2.8  Subroutine  for  Extraction  of  Element  Displacements 
from. the  Global  Displacement  Vector  (XTRACT) 

CALL  XTRACT (LNUM,  NDE,  ELQ) 

LNUM  -  A  positive  scalar  integer  equal  to  user’s  global 
element  number 
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NDE  -  A  scalar  integer  equal  to  the  total  number  of  degrees 
of  freedom  in  the  element 

ELQ  -  A  floating  point  vector  at  least  NDE  words  long. 

The  entries  of  ELQ  are  undefined  when  XTRACT  is 
called. 

If  the  analyst  desires  to  calculate  stress  or  strain  dis¬ 
tributions  in  his  structure,  he  commonly  uses  transformations 
between  stress  or  strain  and  nodal  displacements: 

£  c  B uel  <r  -  Ee 

where  £*,  O'  are  respectively  vectors  of  strain  and  stress  com¬ 
ponents  at  selected  points  in  the  element  domain,  and  where 
Me|  is  the  vector  of  element  nodal  displacements,  i.e.,  a  subset 
of  the  global  vector  U  .  Subroutine  XTRACT  selects  the  correct 
subset  Ue  I  out  of  U  ,  based  upon  the  information  contained 
in  the  user's  master  assembly  list.  The  values  of  lL€f  are 
placed  in  the  argument  vector  ELQ. 
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SECTION  i| 


A  SAMPLE  FEABL  PROGRAM 

In  order  to  illustrate  further  the  data  location  and 
process  sequence  interfacing,  sample  user  programs  will  be  de¬ 
veloped  for  analysis  of  the  truss  structure  shown  in  Figure  12. 
The  structure  consists  of  1 6  bars  and  18  degrees  of  freedom  in 
the  XY  plane. 


FIGURE  12 


Local  number  conventions  for  the  (n) th  typical  bar  element  are 
shown  at  the  right.  With  these  conventions,  the  stiffness  matrix 
of  the  typical  bar  element  is  given  by: 
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(  Sy  ntme^Tl  c) 
c  os3oc 

sine*  cos  oC.  sinloC^ 

where: 

'  =  ■/  (XZ-X1)Z+  (Ya.-Vi)2*' 
cas  ^  ) 

s<*/7  ^  =  ■—  ("  V*  -  Yj ) 

and  where  E,  A  are  the  bar’s  modulus  and  cross  section  area. 

Also,  once  the  displacements  u^,  u^,  u^,  u^  are  known  for  the 
element,  its  elongation  may  be  calculated  as: 

<f  c  (u3~U1)  cosed  +•  (u4  ~Ui  )  Strt* 

and  the  load  in  the  bar  is  then  given  by: 

P  =  EAS/j? 

where  and  oC  are  defined  as  above. 

The  user  decides  to  read  the  properties  and  nodal  coordi¬ 
nates  for  each  element,  each  time  they  are  required  (a  rather 
inefficient  procedure).  Reading  is  to  be  done  by  the  generator 
subroutine,  rather  than  in  MAIN;  however,  the  card  reader  device 
code  will  be  established  in  MAIN.  Therefore,  the  user  programs 
his  stiffness  matrix  and  stress  generator  subroutines  as  follows: 


i 


COS  oi 
sin d(  cot  oi 
-  c OSZ0( 

-sindcotoi 


sinx  oi 


-  stnxeosoi 
-  sin** 
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Stiffness  Matrix  Generator 


SUBROUTINE  BARK  (ELK,  ELQ) 

DIMENSION  ELK  (4,  4),  ELQ  (4) 

COMMON  /I0/  KR,  KW,  KP,  KT1,  KT2 ,  KT3 
C  INTERFACE  WITH  10  CONTROL  PARAMETERS  IS  OPTIONAL 
91  FORMAT  (6E10.3) 

READ  (KR,  91)  XI,  Yl,  X2,  Y2,  E,  A 
C  CALCULATE  BAR  LENGTH 

BARL=SQRT  ((X2-X1)**  2+(Y2-Yl)*»  2) 

C  CALCULATE  SINE  AND  COSINE 
S=( Y2-Y1) /BARL 
C=(X2-X1')  /BARL 

C  CALCULATE  ENTRIES  IN  LOWER  TRIANGLE  OF  ELK-ASEMBL  DOES 
C  NOT  USE  UFPER  TRIANGLE 

ELK  (1,  1) =E*A*C*C/BARL 
ELK  (2,  1)=E*A*S*C/BARL 

etc. 

ELK  (4,  4)=E#A*S*S/BARL 

C  ESTABLISH  ZERO  ELEMENT  EQUIVALENT  N0DAL  FORCES 
DO  10  1=1,  4 
10  ELQ  (I)=0 . 

RETURN 

END 

Stress  Generator 

SUBROUTINE  BARF  (LNUM,  ELQ) 

DIMENSION  ELQ  (4) 

COMMON  /I0/  KR,  KW,  KP,  KT1,  KT2 ,  KT3 

91  FORMAT  (6E10.3) 

92  FORMAT  ( 21H0BAR  FORCE  IN  BAR  N0. ,  14,  2H  =,  E10.3,  3H  LB) 
READ  (KR,  91)  XI,  Yl,  X2,  Y2,  E,  A 

BARL  =  SQRT  ((X2-X1)#*  2+( Y2-Y1)**2) 

S=( Y2-Y1) /BARL 
C=(X2-X1)/BARL 

F0RCE=E*A*(C#(ELQ(3)  -  ELQ  (1))  +  S*  (ELQ(4)-ELQ  (2)))/BARL 

WRITE  (KW,  92)  LNUM,  FORCE 

RETURN 

END 

Subroutines  BARK  and  BARF  are  ready  to  use  in  conjunction  with 
the  FEABL  software.  The  user  now  begins  the  construction  of 
his  MAIN  program,  one  stage  at  a  time. 
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Stage  1:  Program  Heading  (Data  Location  Interface) ,  Device 
Code  and  Problem  Size  Establishment 


The  user  estimates  that  a  1000-word  /DATA/  vector  will 
be  more  than  adequate  for  the  problem.  Required  device  codes 
are  the  card  reader  (5)  and  printer  (6). 

C  MAIN  PR0GRAM  F0R  S0LUTI0N  0F  TRUSS  PR0BLEM 
DIMENSI0N  REAL(IOOO),  INTGR(IOOO) 

DIMENS I0N  ELK(4,  4),  ELQ(4) 

C0MM0N  /I0/  KR,  KW,  KP,  KT1,  KT2 ,  KT3 
C0MM0N  /SIZE/  NET,  NDT 

C0MM0N  /BEGIN/  IC0N,  IK0UNT,  ILNZ,  IMASTR,  IQ,  IK 
C0MM0N  /END/  LC0N ,  LK0UNT,  LLNZ,  LMASTR ,  LQ,  LK 
EQUIVALENCE  (REAL(l),  INTGR(l) ) 

KR=5 

KW=6 

NET=l6 

NDT=l8 

CALL  SETUP  (1000,  5,  80) 

C  END  0F  STAGE  1 

The  user  has  called  for  space  for  five  constraints:  master 
displacements  2  (after  a  -45°  rotation  of  1  and  2),  7,  8,  13 
and  14.  The  assembly  list  must  allow  room  for  16  elements  x 
(4  DOF  plus  1  pointer  per  element)  =  80  words.  The  first 
declaration: 

DIMENSI0N  REAL(IOOO),  INTGR(IOOO) 
is  also  duplicated  and  placed  in  each  FEABL  subroutine. 

Stage  2:  Assembly  List  Input  and  Organization  of  K 

The  user  recognizes  that  the  pointers  will  occupy  locations 
IMASTR  to  IMASTR+15  (see  Section  2,2.4).  He  chooses  to  write 
specific  assignment  instructions  for  each  element: 

INTGR( IMASTR) =IMASTR+l6  (Pointer  for  1st  element) 

INTGR( IMASTR+16) =1 

INTGR( IMASTR+17 ) =2 

INTGR( IMASTR+18) =7 

INTGR ( IMASTR+19 ) =8 

• 
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Stage  1:  Program  Heading  (Data  Location  Interface) ,  Device 
Code  and  Problem  Size  Establishment 


The  user  estimates  that  a  1000-word  /DATA/  vector  will 

be  more  than  adequate  for  the  problem.  Required  device  codes 

are  the  card  reader  (5)  and  printer  (6). 

C  MAIN  PR0C-RAM  F0R  S0LUTI0H  0F  TRUSS  PR0BLEM 
DIMENSI0N  REAL(IOOO),  INTGR(IOOO) 

DIMENSI0N  ELK(4,  4),  ELQ( 4 ) 

C0MM0N  /I0/  KR,  KW,  KP,  KT1,  KT2,  KT3 
C0MM0N  /SIZE/  NET,  NDT 

C0MM0N  /BEGIN/  IC0N,  IK0UNT,  ILNZ,  IMASTR,  IQ,  IK 
C0MM0N  /END/  LC0N,  LK0UNT,  LLNZ,  LMASTR ,  LQ,  LK 
EQUIVALENCE  (REAL(l) ,  INTGR(l)) 

KR=5 

KW-6 

NET=l6 

NDT=l8 

CALL  SETUP  (1000,  5,  80) 

C  END  0F  STAGE  1 

The  user  has  called  for  space  for  five  constraints:  master 
displacements  2  (after  a  -45°  rotation  of  1  and  2),  7,  8,  13 
and  14.  The  assembly  list  must  allow  room  for  16  elements  x 
(4  DOF  plus  1  pointer  per  element)  =  80  words.  The  first 
declaration: 

DIMENSI0N  REAL(IOOO),  INTGR(IOOO) 

is  also  duplicated  and  placed  in  each  FEABL  subroutine. 

Stage  2:  Assembly  List  Input  and  Organization  of  K 

The  user  recognizes  that  the  pointers  will  occupy  locations 

IMASTR  to  IMASTR+15  (see  Section  2.2.4).  He  chooses  to  write 

specific  assignment  instructions  for  each  element: 

INTGR( IMASTR) =IMASTR+l6  (Pointer  for  1st  element) 

INTGR( IMASTR+16 ) =1 
INTGR( IMASTR+17 ) =2 
INTGR ( IMASTR+ 1 8 )  =  7 
INTGR ( IMASTR+19 ) =8 


42 


INTGR( IMASTR+7 ) =IMASTR+44  (pointer  for  8th  element) 

INTGR ( IMASTR+4  4  )  =7 
INTGR  ( IMASTR+4  5 )  =8 
INTGR ( IMASTR+4 6 ) =9 

INTGR  ( IMASTR+4  7 )  =10 

« 

• 

INTGR ( IMASTR+15) =IMASTR+76  (pointer  for  16 th  element) 

INTGR ( IMASTR+7  6 )  =  7 
INTGR (IMASTR+7 7) =8 
INTGR (IMASTR+7 8) =15 

INTGR(IMASTR+79)=l6  (80th  location  in  master  assembly  list) 

At  this  point,  it  is  advisable  to  dump  the  assembly  list  for 
debugging  purposes  if  there  is  any  doubt  about  its  accuracy. 

Dumping  may  be  done  by: 

WRITE  (KW,  5)  (INTGR(I),  I=IMASTR,  LMASTR) 

5  F0RMAT  (IX,  10110) 

Finally: 

CALL  0RK( 1000) 

C  END  0F  STAGE  2 

Stage  3:  Generation  and  Assembly  of  Element  Properties 

In  this  stage,  the  user  merely  invokes  the  appropriate 

subroutines  in  a  loop. 

D0  10  LNUM=1 ,  NET 
CALL  BARK  (ELK,  ELQ) 

CALL  ASEMBL  (LNUM,  4,  ELK,  ELQ) 

10  C0NTINUE 

C  END  0F  STAGE  3 

Stage  4:  Rotation  Transformations 

In  the  present  problem,  only  the  lower  left  node  (Figure  12) 

requires  rotation.  Master  degrees  of  freedom  number  1  and  2 

must  be  rotated  by  -45  degrees. 

CALL  R0TATE(1,  .1,  2,  0,  -45.,  0.,  0.) 

t  ft 

Quantities  for  3-D  problems  are  not  used 

C  END  0F  STAGE  4 
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Stage  5:  Boundary  Conditions 

Displacements  2,  7,  8,  13,  1 4  are  constrained,  and  the 

constraint  vector  starts  in  INTGR(l)  (IC0N=1).  Therefore: 

INTGR(l) =2 
IHTGR( 2)=7 
INTGR(3)=8 
INTGR(4)=13 
INTGR(5)=14 

The  constrained  displacements  are  all  prescribed  to  be  zero. 

Only  the  nonzero  prescribed  forces  need  be  input;  these  are  1,000 

lb.  each  at  the  15th  and  17th  degrees  of  freedom.  Therefore: 

C  PRESCRIBED  DISPLACEMENTS 
ELAL(IQ+l)-0 . 

REAL(IQ+6)=0. 

REAL(IQ+7)=0. 

REAL(IQ+12)=0. 

REAL(IQ+13)=0. 

C  ACCUMULATE  PRESCRIBED  F0RCES 
REAL( IQ+14 ) =REAL( IQ+14 ) +1000 . 

REAL( IQ+16 ) =REAL( IQ+16 )+ 1000 . 

CALL  BC0N 
C  END  0F  STAGE  5 

Stage  6:  Choleski  Solution 

CALL  FACTPD 

CALL  SIMULQ( ENERGY) 

C  END  0F  STAGE  6 

Stage  7:  Inverse  Rotation  to  Obtain  Global  Displacements  in 
Global  Cartesian  Coordinates 
A  rotation  of  -45°  was  performed  at  node  1  before  the 
boundary  conditions  were  imposed.  Since  all  element  calculations 
are  done  with  respect  to  the  unrotated  XY  axis  system,  this  rota¬ 
tion  must  be  reversed  before  element  stresses  are  calculated: 

CALL  R0TATE(-1,  1,  2,  0,  45.,  0.,  0.) 

C  END  0F  STAGE  7 

Stage  8:  Calculation  of  Element  Stresses 

In  this  case,  bar  forces  are  to  be  calculated.  Again, 
the  user  merely  invokes  the  proper  subroutines: 
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D0  20  LNUM=1,  NET 

CALL  XTRACT  (LNUM,  4,  ELQ) 

CALL  BARF  (LNUM,  ELQ) 

20  C0NTINUE 
C  END  0F  PR0G 

ST0P 
END 

If  the  reader  has  grasped  the  material  presented  up  to 
this  point,  he  now  has  enough  familiarity  with  FEABL  to  use  it, 
albeit  somewhat  inefficiently.  However,  some  additional  degree 
of  sophistication  is  desirable  for  the  instructions  which  input 
problem  data  such  as  the  assembly  list.  It  is  apparent  that 
a  straightforward  net  of  instructions  such  as  that  given  in  the 
sample  program  is  quite  cumbersome,  especially  for  problems  in¬ 
volving  more  elements  with  more  degrees  of  freedom  per  element. 
Improved  programming  techniques  are  discussed  in  the  next 
section. 
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SECTION  5 


ADVANCED  TECHNIQUES  WITH  PEABL 

5.1  Efficient  Programming;  for  Large  Problems 

When  a  structure  is  to  be  analyzed  with  a  large  number 
of  elements  and  many  degrees  of  freedom,  there  often  occurs  a 
definite  trend  toward  regularity  in  the  element  set.  The  analyst 
should  recognize  two  distinct  forms  of  regularity,  and  he  should 
be  prepared  to  take  advantage  of  each  in  writing  his  MAIN  program. 

First,  consider  the  truss  structure  shown  in  Figure  13, 
with  all  bays  having  the  same  dimension  vertically  and  horizontally. 


FIGURE  13 


US 


The  bar  elements  may  be  said  to  possess  a  degree  of  regularity, 
in  that  there  are  only  four  different  types  of  elements  in  the 
structure: 

1.  Horizontal  bars 

2.  Vertical  bars. 

3.  145°  diagonals. 

4.  -45°  diagonals. 

Since  an  element  stiffness  matrix  depends  upon  nodal  coordinates 
only  through  differences  in  the  coordinate  values  (see  beginning 
of  Section  4),  there  will  occur  only  four  independent  element 
stiffness  matrices  for  the  structure  of  Figure  13:  one  for  each 
of  the  element  types  listed  above.  Thus,  for  example,  the  element 
stiffness  matrix  for  a  typical  horizontal  bar  may  be  generated 
once  and  assembled  repeatedly  (32  times  for  the  present  case). 

The  technique  may  be  repeated  again  for  the  30  vertical  bars, 

24  +45°  diagonals  and  24  -45°  diagonals,  achieving  considerable 
savings  in  execution  time. 

The  second  form  of  regularity  involves  the  way  in  which 
displacements  and  elements  may  be  numbered.  Often  an  element 
set  is  topologically  equivalent  to  a  rectangular  or  square  array, 
even  if  it  is  not  geometrically  regular.  For  example,  suppose 
the  trapezoidal  continuum  shown  in  Figure  14(a)  is  to  be  analyzed 
in  plane  stress  using  32  3-node  triangle  elements.  For  numbering 
purposes  the  element  and  displacement  sets  are  topologically 
equivalent  to  the  square  net  shown  in  Figure  14(b).  It  is  then 
possible  to  generate  the  assembly  list  for  the  structure  by  a 
double  D0  loop.  Each  element  .is  ’'located”  in  the  structure  via 
an  intersection  of  one  of  the  element  strings  LX  and  one  of  the 
element  strings  LY.  The  element  number  and  its  master  displace¬ 
ment  numbers  may  be  generated  from  the  values  of  LX  and  LY. 

Taking  advantage  of  numbering  regularity  enables  the  user  to 
input  large  sections  of  the  assembly  list  with  relatively  few 
F0RTRAN  instructions. 


FIGURE  14 


(a.) 


(W 


A  Numbering  Strategy  for  Planar  Problems 


Suppose  the  set  of  MxN  rectangular  plane  stress  el  aments 
shown  In  Figure  15  is  to  be  assembled,  where  M,  N  may  be  quite 
large.  If  the  numbering  strategy  shown  in  the  figure  is  adopted, 
beginning  with  element  1  and  displacements  1,  2  at  the  lower 
left  corner  and  ending  with  element  MN  and  displacements 
2(M+1)  (N+l)-l,  2(M+1)  (N+l)  at  the  upper  right  corner,  then 
each  element  number  can  be  given  as  a  function  of  its  LY,  LX 
string  coordinates; 

LNUM=LY+N*(LX-1)  where  1<LY<N  and  1<LX<M 
In  order  to  create  similar  functions  for  the  element  assembly 
list,  a  local  numbering  convention  must  be  adopted.  If  the 
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FIGURE  15 

convention  shown  is  Figure  15  is  used,  the  global  DOF  numbers 

may  be  calculated  as  follows:  Local  displacements  1,  2  can  be 

considered  to  lie  on  strings  LX,  LY;  hence: 

(Local  #1)*2*LY-1+2«(N+1)«(LX-1) 

(Local  #2)=(Local  #1)+1 

Local  displacements  7,  8  follow  directly  from  the  master  numbering 
scheme: 

(Local  #7)=(Local  #2)+l=(Local  #l)+2 
(Local  #8)=(Local  #7)+l=(Local  #1)+ 3 

Local  displacements  3,  4  are  shifted  one  LX  string  to  the  right; 

hence  they  are  2#(N+1)  ahead  of  1,  2: 

(Local  #3)*(Local  #1)+2*(N+1) 

(Local  #4)=(Local  #l)+2*(N+l)+l=(Local  #2)+2*(N+l) 
and  finally: 

(Local  #5)=(Local  n)+2*(N+l)+2=(Local  #7)+2*(N+l) 

(Local  #6)=(Local  #l)+2*(N+l)+3=(Local  #8)+2*(N+l) 

With  the  algorithms  derived  above,  it  is  quite  easy  to 

develop  an  efficient  procedure  to  input  the  entire  assembly  list 

into  the  /DATA/  vector.  The  reader  will  recall  that  the  pointer 
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for  element  number  LNUM  is  stored  in  INTGR( IMASTR+LNUM-1) . 

One  additional  Integer  variable,  called  NEXT,  is  required.  NEXT 
is  to  be  incremented  after  completion  of  the  input  for  one  ele¬ 
ment,  so  that  the  next  value  of  NEXT  is  the  address  of  the  next 
available  location  in  the  master  assembly  list,  i.e.,  NEXT= 
the  value  oi  the  pointer  for  the  next  element.  The  entire  assembly 
list  is  then  generated  and  input  by  the  following  set  of  P0RTRAN 
instructions : 

(Values  of  M  and  N  defined  elsewhere  in  program) 

C  INITIALIZE  P0INTER  VALUE 
NEXT=IMASTR+NET 

C  L00P  0VER  ELEMENT  STRINGS 
D0  20  LX-1,  M 
D0  20  LY=1,  N 
LNUM=LY+N* (LX-1) 

I PTR= IMASTR+LNUM-1 

C  ESTABLISH  P0INTER  F0R  THE  ELEMENT 
INTGR( IPTR)=NEXT 

C  CALCULATE  JD1-1  LESS  THAN  MASTER  N0.  0F  1ST  D0F,  JD3=1 

C  LESS  THAN  MASTER  N0.  0F  3RD  D0F 
JD1=2*(LY-1)+2*(N+1)*(LX-1) 

JD3=JD1+2*(N+1) 

C  ASSIGN  DOF  MASTER  NOS.  TO  LOCATIONS  NEXT,  NEXT  +  1,...,  NEXT  +  7 

C  IN  /DATA/  VECT0R 
J=0 
K=0 

D0  10  1=1,8 
INDEX=NEXT+I-1 

IF( I  .GT.  2  .AND.  I  .LT.  7)  G0  T0  5 

J=J+1 

JD=JD1 

II=J 


The  last  IP  statement  above  is  a  good  form  of  insurance  for  cases  f 

in  which  the  user  may  have  overestimated  the  number  of  words  j 

required  for  his  master  assembly  list.  1 

] 

A  similar  scheme  for  structures  consisting  of  triangle  > 

elements  (Figure  1*0  can  be  developed  if  the  "right-side-up"  1 

and  "upside-down"  elements  are  treated  as  separate  sets. 

The  automatic  generation  technique  derived  above  was  de-  j 

veloped  for  a  structure  with  complete  topological  regularity; 
however,  it  may  be  extended  to  cover  large  portions  of  less  regular  ' 

structures  with  relatively  little  additional  programming.  There 
is  also  a  hidden  advantage:  the  numbering  strategy  adopted  in  \ 

Figure  15  is  not  only  easy  to  produce,  but  a3.so  minimizes  the  f 

population  of  the  master  stiffness  matrix  if  none  of  the  elements  ■ 

have  mid-side  nodes.  j 

i 

5.3  A  Short  Note  on  Automatic  Data  Generation  \ 

~~  J,ni  t 

% 

The  reader  will  recall  that  the  hypothetical  user  in  J 

Section  *1  chose  to  read  in  his  element  coordinates  and  basic  1 

properties  from  data  cards.  This  procedure  can  be  time-consuming 
and  expensive  for  problems  involving  large  numbers  of  elements, 
especially  when  the  element  set  geometry  is  not  regular  in  the 
sense  of  Subsection  5.1  and  Figure  13.  However,  it  often  happens 
that  the  element  set  geometry  is  regular  in  the  sense  that  the 

element  coordinates  may  be  calculated  from  a  general  algorithm  i 

based  upon  the  element  string  concept  discussed  in  the  previous 

subsection,  i.e.: 

D0  20  LX=1,  M 
D0  20  LY=1,  M 
LJUM=LY+N#(LX-1) 

(M(LX,  LY)) 

• 

• 

where  X={X1,  Yl,  X2,  Y2,...}  is  the  element  coordinate  vector 
and  f  is  a  floating  point  function  of  LX  and  LY.  An  algorithm 
of  this  type  may  be  used  to  generate  the  data  base,  as  it  is 
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needed,  for  calculation  of  element  stiffness  matrices  and  element 
stresses.  This  results  in  the  trade-off  of  a  slight  increase 
in  execution  time  and  some  decrease  in  required  storage  space, 
since  it  is  not  necessary  to  carry  large  vectors  of  global  nodal 
coordinates  in  core. 

5.4  A  Handy  Trick 

There  often  occur  plane  elasticity  problems  in  which  a 
large  number  of  degrees  of  freedom  are  to  have  prescribed  displace¬ 
ments.  For  example,  suppose  a  rectangular  domain  such  as  the 
one  shown  in  Figure  15  is  to  be  analyzed,  and  that  all  four  edges 
of  the  domain  are  clamped.  This  means  that  25  to  30  per  cent 
to  the  total  degrees  of  freedom  will  have  prescribed  displacements. 

A  significant  amount  of  execution  time  may  be  saved  by 
modifying  the  numbering  scheme  shown  in  Figure  15,  so  that  the  edge 
degrees  of  freedom  have  the  largest  global  numbers.  Let  NDT 
be  the  total  number  of  degrees  and  NFT  be  the  total  number  of 
unconstrained  degrees.  Then  the  modified  number  scheme  assigns: 

1,  2,...,  NFT 

to  the  unconstrained  degrees  and: 

NFT+1,  NFT+2 , . . . ,  NDT 

to  the  degrees  along  the  edges  of  the  structure.  Program  stages 
1  through  5  are  completed  in  standard  fashion. 

However,  just  before  factoring  K  ,  the  user  may  fool 

FEABL  by  inserting: 

ITEMP=NDT 

NDT=NFT 

in  his  MAIN  subprogram.  Subroutines  FACTPD/FACTSD  and  SIMULQ 
will  then  solve  only  for  the  unknown  displacements  1,  2,...,  NFT. 

At  the  beginning  of  Stage  7,  when  the  full  displacement  vector 
may  be  required  again,  the  user  inserts: 

NDT=ITEMP 
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The  real  value  of  this  trick  depends  on  a  trade-off  between 
execution  time  and  core  storage.  The  modified  numbering  scheme 
will  result  in  a  requirement  for  additional  storage  space  for  K  , 
over  what  is  needed  by  the  numbering  scheme  discussed  in  Subsection  5.2. 
The  excess  requirement  is  given  by: 

AS»(NDT-NFT)  (NDT-B) 

where  B  is  the  average  semi-bandwidth  of  the  unconstrained  Dart 

of  K  . 


5.5  Time-Saving  Techniques  for  Design  Studies 

The  "design  study"  approach  to  finite  element  analysis 
may  take  any  of  the  following  forms: 

1.  Consideration  of  a  number  of  loading  environments 
applied  to  a  unique  structure  with  unique  displacement 
boundary  conditions. 

2.  Analysis  of  a  unique  structure  under  various  environments 
in  which  prescribed  displacements  are  changed,  as 

well  as  prescribed  forces. 

3.  Consideration  of  variations  in  the  structure  itself 
to  meet  a  given  environment. 

The  first  two  categories  are  self-explanatory.  An  example  of 
the  third  might  be  an  analysis  of  a  truss  structure  in  which 
the  compression  bars  are  checked  for  buckling: 

P  <  P=  7 rzGl/jz 

If  any  bar  were  found  to  be  under  a  compressive  load  greater 
than  its  P  ,  the  program  might  incorporate  an  algorithm  for 

w  X 

redesigning  the  bar  and  re-analyzing  the  new  structure. 

The  above  categories  of  problems  can  be  studied  using 
FEABL  with  only  minor  modifications.  No  changes  are  required 
in  the  data  location  interfaces  and  if  a  change  to  the  process 
sequence  interface  is  required,  it  involves  only  an  intelligent 
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application  of  the  rules  presented  above.  However,  the  analyst 
must  be  familiar  with  programming  techniques  required  for  com¬ 
munication  between  his  computer's  core  and  external  storage  devices 
(system  disks,  drums,  tape  units).  These  techniques  vary  from 
installation  to  installation  and  will  not  be  discussed  in  specific 
terms  here.  The  analyst  should  consult  the  operating  manuals 
applicable  to  the  system  which  he  will  be  using. 

Load  environment  case  study  is  the  easiest  type  of  multi¬ 
solution  problem  to  handle.  Since  the  displacement  boundary 
conditions  are  not  varied,  only  one  assembly  and  factoring  of 
the  master  stiffness  matrix  need  be  done.  The  factored  form 
of  K  is  held  in  the  /DATA/  vector  while  each  prescribed  vector 
is  formed  and  the  displacements  are  solved  for.  The  quantity: 

«> 

A 

Q =  Assembled  Element  Equivalent  Nodal  Forces 

£  A 

and  the  prescribed  displacements  IL  D  must  be  saved.  These  two 
quantities  are  found  in  the  Force/Displacement  Vector: 


in  the  /DATA/  vector  just  after  BC0N  is  called.  (The  prescribed 
external  loads  Q  F,  on  the  other  hand,  are  not  really  needed 
until  SIMULQ  is  called.)  Figure  16  illustrates  the  modified 
FEABL  process  sequence  whi^h  will  accomplish  this  result. 

The  second  category  involves  essentially  the  same  techniques 
as  the  first.  However,  in  this  case  the  data  which  must  be  saved 
in  external  storage  are: 

1.  Cl  E  -  Assembled  element  equivalent  nodal  forces. 

2.  K  -  Assembled  master  stiffness  matrix  prior  to 

boundary  condition  application. 
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FIGURE  16 

For  each  case  the  user  must  input  the  contents  of  the  constraint 

A  A  ^ 

vector,  fetch  Q  £,  accumulate  and  reset  fetch  K 

and  finally  re-enter  the  standard  FEABL  process  at  the  point 
where  BC0N  is  called. 

Significant  time  savings  are  gained  in  the  third  problem 
category  if  the  design  changes  involve  only  a  few  elements  in 
a  structure  composed  of  a  large  number  of  elements.  Suppose 
that  the  displacement  solution  has  been  found  from  the  f'-rce- 
displacement  relations  for  the  initial  structure: 

K  u  «  Qf*Q* 

During  calculation  of  the  element  stresses,  it  is  discovered 
that  design  changes  are  required  in  one  or  more  elements;  these 
changes  will  appear  in  the  force-displacement  relations  as  a 
modification  aK^  of  the  master  stiffness  matrix  and  (possibly) 
a  modification  A  Q  E  of  the  assembled  element  nodal  force  vector. 
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The  new  force-displacement  relations: 

(K  ♦  aK)  «'  =  Qr  <•  Q*  t  LQ* 

must  then  be  solved  for  the  modified  displacements  vf  . 

To  program  the  design  change  technique  with  FEABL,  the 
user  must  take  the  following^  nonstandard  actions: 

1.  Output  K  and  <?P+<?  ^  to  external  storage  after 
assembly  and  input  of  prescribed  quantities,  but  before 
calling  BC0N. 

2.  After  the  initial  displacement  solution  U  has  been 
obtained,  transfer  it  to  temporary  core  storage  out¬ 
side  the  /DATA/  vector.  (The  Force/Displacement  Vector 
area  in  the  /DATA/  vector  will  be  required  for  accu- 

A 

mulation  of  A$g.)  Zero  the  Force/Displacement  Vector 
and  the  Master  Stiffness  Matrix  Array. 

3.  Calculate  stresses  from  V-  ,  element-by-element.  (The 
user  will  need  his  own  version  of  XTRACT  to  extract 
the  proper  #el  from  temporary  core  storage.)  When^ 

a  design  change  is  required,  calculate  aUs1  and  A#el 
and  accumulate  them  using  FEA3L  subroutine  ASEMBL. 

4.  Apply  rotation  transformations  (if  any)  to  aK  and 

A  A 

A reset  entries  of  A£?g  to  zero  where  displacements 
are  prescribed. 

5.  Fetch  K  and  from  external  storage  and  accumu¬ 

late  them  to  A  K  and  AjSf^.  Apply  boundary  conditions 
and  solve  for  vf  . 

The  above  process  may  be  repeated  in  an  iterative  design 
process,  with  K  +A  K,  Ct  +A$*  at  the  beginning  of  each 

new  design  step  playing  the  role  of  the  initial  values  K  >$p+$  • 
This  procedure  is  referred  to  as  iterative  updating,  and  is  also 
applicable  to  nonlinear  elastic  and  elastic-plastic  analysis 
of  continua.  In  these  types  of  analysis  the  "design  change" 
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results  from  following  a  material  stress-strain  curve  and/or 
testing  the  satisfaction  of  a  yield  conditon  (e.g,,  the  Mises- 
Hencky  criterion). 


Substructurim 


with  FEABL 


The  FEABL  software  system  has  been  designed  primarily 
for  in-core  solutions.  Up  to  1,500  degrees  of  freedom  can  be 
handled  on  currently  available  hardware  with  500  KBYTE  (125  KVJord) 
memory.  FEABL* s  in-core  capability  may  be  extended  by  means 
of  the  substructuring  technique  outlined  briefly  here. 

Figure  17  illustrates  a  domain  which  has  been  divided 
into  a  small  number  of  substructures.  The  desired  stress  solution 
accuracy  is  on  a  scale  much  smaller  than  the  substructure  dimension, 


Sub  struct 

#i 


Sub  struct  Substriurt 


-»eSircov. 

accuracy 


Substructure  divided 
into  ordinary  elements 


FIGURE  17 

so  each  substructure  is  subdivided  further  into  ordinary  elements. 
Let  the  subscripts  I  and  B  refer,  respectively,  to  the  interior 
and  boundary  degrees  of  freedom  in  a  substructure.  Then  the 


force-displacement  relations 
into:  1^- 

K 


for  the  substructure  may  be  partitioned 


The  interior  degrees  of  freedom  may  be  eliminated  by  the  process 
of  static  condensation,  which  transforms  the  substructure  force- 
displacement  relations  to  the  form: 

~  (  ^Bd  “  ^bi  ^ix  ^  ~  “  ^Br^rr  =  Qt 

Each  substructure  may  now  be  considered  as  a  "superelement" 
having  degrees  of  freedom  only  on  the  interelement  boundaries, 
as  shown  in  Figure  18.  The  quantities  and  are  then 


FIGURE  18 

assembled  in  the  same  manner  as  the  stiffness  matrices  and  equiva¬ 
lent  nodal  force  vectors  for  ordinary  elements.  After  the  displace¬ 
ment  solution  on  the  substructure  boundaries  has  been  obtained, 
the  original  force-displacement  relations  for  each  substructure 
may  be  used,  with  U  ^  as  prescribed  displacements,  to  obtain 
the  interior  solution. 
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In  addition  to  extension  of  FEABL’s  problem  size  capability, 
the  substructuring  technique  may  also  be  used  to  reduce  roundoff 
error  and  to  save  execution  time.  All  of  these  improvements 
depend  upon  optimization  of  the  relative  numbers  of  substructures 
and  of  ordinary  elements  within  the  substructures.  Use  of  the 
substructuring  technique  with  FEABL  software  will  be  discussed 
in  detail  in  a  future  publication. 

5.7  Error  Estimation  Methods 

Rigorous  mathematical  proofs  exist  showing  that  a  finite 
element  analysis  for  the  stress  and  displacement  distributions 
in  a  structure  converges  to  the  exact  solution  as  the  number  of 
elements  in  a  given  domain  is  increased,  provided  only  that  certain 
easily  satisfied  restrictions  are  obeyed  (Ref. 3).  However,  since 
all  calculations  done  in  a  digital  computer  are  imprecise,  errors 
due  to  roundoff  will  occur.  It  has  been  shown  (Ref.  l\ )  that 
roundoff  error  increases  in  an  extremely  complex  way  as  the  number 
of  elements  is  increased.  Therefore,  it  is  advisable  to  resort 
to  some  numerical  method  which  will  produce  a  reasonable  estimate 
of  the  error,  and  which  does  not  depend  upon  the  details  of  what 
types  of  elements  are  used  or  what  boundary  conditions  are  applied 
to  the  structure.  Four  methods  are  presented  here. 

5.7«1  Irons*  Energy  Variance  Criterion 

The  energy  criterion  proposed  by  Irons  (Ref.  5)  is  the 
most  economical,  requiring  no  more  time  than  the  calculation 
of  the  rounding  error  parameter  discussed  in  Subsection  3.2.6. 

The  diagonal  entries  of  the  master  stiffness  matrix  must 
be  saved,  either  in  temporary  core  storage  or  on  an  external 
unit.  After  the  displacement  solution  has  been  obtained,  the 
energy  variance  can  be  calculated  from: 
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Z  *v  B  not  Jo/ 


z 

1=1  J 


SK‘i“'“4 


where  p  is  the  computer  precision  in  decimal  places  and  B  is 

the  average  semi-bandwidth  of  K  .  The  value  of  B  can  be  obtained 

easily  from  the  address  index  parameters*: 

B=FL0AT( LK+l-IK ) /FL0AT ( NDT ) 

One  half  the  value  of  the  denominator  of  the  energy  variance 
expression  is  provided  to  the  user  by  FEABL  subroutine  SIMULQ 
in  the  argument  ENERGY  of  that  subroutine. 

Some  care  must  be  exercised  in  applying  Irons  ’  energy 
criterion.  The  user  will  obtain  unrealistically  large^alues 
of  (f  in  problems  in  which  the  structure  is  constrained  very 
lightly,  and  in  which  large  amounts  of  "self-energy"  can  be  stored 
when  a  single  degree  of  freedom  is  displaced.  The  cantilever 
beam  is  a  good  example  of  a  structure  to  which  the  Irons’  criterion 
cannot  be  applied. 

5.7.2  The  Residual  Force  Method 

Calculation  of  residual  forces  provides  a  more  detailed 
picture  of  the  distribution  of  errors  through  the  structure. 

The  master  stiffness  matrix  may  be  saved  in  external  storage  for 
this  purpose  Immediately  after  calling  BC0N.  Approximate  values 
of  the  reaction  forces  at  degrees  where  displacements  were  pre¬ 
scribed  may  be  obtained  as  well  by  saving  K  just  prior  to  calling 
BC0N. 


FL0AT  is  an  IBM  library  function  which  converts  integers  to 
floating  point  numbers. 
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Let  U  p  be  the  approximate  displacement  solution  obtained 
from  the  constrained  force-displacement  relations: 

KFF  0  fQp'KFD^ 

lo  ijlij  l  fi* 


Then: 


is  the  approximate  force  vector.  When  K  has  been  returned 
to  core  after  the  displacement  solution  has  been  obtained,  Q* 
can  be  calculated  by  the  following  algorithm: 

DIMENSI0N  R(500) 

C  VECT0R  R  MUST  ALL0W  EN0UGH  ST0RAGE  P0R  THE  PULL 
C  P0RCE  VECT0R  Q-STAR  IF  REACTI0N  F0RCES  ARE 
C  DESIRED.  0THERWISE,  EN0UGH  P0R  THE  UNC0NSTRAINED 
C  DEGREES  0F  PREED0M  MUST  BE  ALL0WED 
C  ALLOCATE  EXTERNAL  PILES 

(Standard  FEABL  process  sequence,  except  where  noted) 


C  K  MATRIX  T0  EXTERNAL  ST0RAGE 
WRITE  (...)  (REAL(I),  I=IK,LK) 

C  P0RCE/DISPL  MUST  ALS0  BE  SAVED  F0R  C0MPARIS0N  LATER 
WRITE  (...)  (REAL(I) ,  I*IQ,  LQ) 

CALL  BC0N 
CALL  FACTPD 

CALL  SIMULQ(STRE) 

• 

(Stress  solution,  if  desired) 

C  RETURN  K  MATRIX  F0R  CALOULATI0N  0F  Q-STAR 
READ  (...)  (REAL(I),  I=IK,  uK) 

C  INITIALIZE  C0NSTRAINED  R0W  P0INTER  AT  1ST  N0NZER0  R0W 
D0  50  II=IC0N,  LC0N 
IF(INTGR(II)  .EQ.  0)  G0  T0  50 
NC=II 
G0  T0  60 
50  C0NTINUE 
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C  INITIALIZE  R  VECT0R  P0INTER 
60  NR=1 

C  L00P  0VER  DEGREES  0F  FREED0M,  SKIPPING  C0NSTRAINED  DEGREES 
D0  70  IR0W=1 s  NDT 
IF( INTGR(MC)  ,NE.  IR0W)  G0  T0  71 
C  UPDATE  R0W  P0INTER 
NC=NC+1 
G0  T0  70 

C  F0RM  SUMMATI0N  0F  K(IR0W,  J)*  U(  J)  ,  J=LNZ  T0  NDT 

71  INIT=ILNZ+IR0W-1 
INIT=INTGR( INIT) 

KK=IK0UNT+IR0W-1 
KK=INTGR( KK) 

SUM=0 . 

D0  72  J=INIT,  IR0W 

KADR=KK+J 

JJ=IQ+J-1 

72  SUM=SUM-'-REAL(KADR)  #  REAL(JJ) 

C  REMAINDER  0F  SUM  MUST  TAKE  K  ENTRIES  FR0M  C0L  IR0W 
INIT=IROW4-l 

IF( INIT  .GT.  NDT)  G0  T0 
D0  73  J=INIT,  NDT 

C  MAKE  SURE  R0W  J  HAS  AN  ENTRY  IN  C0L  IR0W 
JJ=ILNZ+J-1 

IF(INTGR(J.J)  .GT.IR0W)  G0  T0  73 
KADR=IK0UNT+J-1 
KADR=INTGR( KADR)+IR0W 
JJ=IQ+J-1 

SUM=SUM+REAL( KADR) *REAL( J J ) 

73  C0NTINUE 

C  PLACE  SUM  IN  NEXT  AVAIL  R  L0CATI0N  AND  UPDATE  R  P0INTER 
71!  R(NR)=SUM 
NR=NR+1 
70  C0NTINUE 

C  RETURN  0RIGINAL  Q  T0  C0RE  F0R  C0MPARIS0N 
READ  (...)  (REAL(I),  I=IQ,  LQ) 

C  CALCULATE  RESIDUAL  F0RCES  Q-(Q-STAR)-Q-R 
C  REINITIALIZE  NC  AND  NR 
NC=II 
NR=1 

D0  80  IR0W=1,  NDT 
IF(INTGR(NC)  .NE.  IR0W)  G0  T0  8l 
NC=NC+1 
G0  T0  80 
81  JJ=IQ+IR0W-1 

R(NR)=REAL( JJ)-R(NR) 

WRITE  (KW,  300)  IR0W,  R(NR) 

500  F0RMAT  ( 20H  RESIDUAL  F0RCE  H0.,l6,  IX,  1H=,  E10.3) 

NR=NR+1 
80  C0NTINUE 
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In  the  above  algorithm,  the  vector  of  residual  forces: 


has  been  calculated.  The  entries  of  Rp  provide  the  detailed  pic¬ 
ture  mentioned  at  the  beginning  of  this  subsection.  Another 
useful  parameter  for  overall  error  measurement  is  the  force  vector 
magnitude  ratio: 


g(R  f) 


Vz. 


Lze<3/-2KB«j): M 


uqpl 
1 9,1 


where  the  summations  extend  only  over  the  unconstrained  degrees 
of  freedom. 

5. 7.- 3  Re-Solution  for  Residual  Displacements 

Although  the  residual  force  vector  is  fairly  easy  to  cal¬ 
culate,  the  interpretation  of  its  meaning  is  not  a  trivial  task. 

The  displacement  solution  must  certainly  be  judged  acceptable 
if,  for  example,  there  are  many  residual  forces  on  the  order 
of  1  lb.  at  degrees  where  forces  of  1,000  lb.  were  applied  originally. 
However,  a  common  situation  in  finite  element  analysis  is  that 
the  applied  force  is  zero  at  many  degrees  of  freedom.  What  does 
a  1  lb.  residual  force  mean  at  these  points?  The  averaged  mea¬ 
sure  rp  presented  in  the  previous  section  relieves  this  detailed 
interpretation  problem  to  some  degree.  However,  the  averaged 
measure  of  greatest  interest  to  the  analyst  is  the  displacement 
vector  magnitude  ratio: 

.  |awf  l 
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Unfortunately,  no  simple  relation  exists  between  dp  and  the  force 
vector  magnitude  ratio  /*p.  If  K  is  an  ill-conditioned  matrix, 
dw  may  in  fact  be  much  larger  than  Y'p. 

If  the  analyst  is  willing  to  spend  some  additional  computing 

time,  the  displacement  residuals  and  an  approximate  calculation 

of  dp  may  be  obtained  by  re-solution.  This  technique  requires 

an  additional  external  storage  file  capable  of  holding  K  .  Just 

before  K  is  returned  to  core  in  the  algorithm  of  Subsection  5.7.2, 

T 

its  factored  form  is  read  into  this  extra  file.  Then, 

picking  up  where  the  previous  algorithm  ended,  the  original  dis¬ 
placement  vector  magnitude  is  calculated,  ldlt  is  returned 
to  core  and  the  contents  of  R  are  transferred  to  the  proper 
locations  in  the  force/displacement  block  of  the  /DATA/  vector. 
Re-solution  is  now  done  simply  by  calling  SIMULQ  again,  after 
which  A U  p  will  be  found  in  the  force/displacement  block.  The 
displacement  vector  magnitude  ratio  obtained  from  this  procedure 
is  actually: 


and  dp*^dp  unless  |  AW  p|  <<  |  U  p*  |  . 

5.7.4  The  Method  of  Rigid  Body  Modes 

A  somewhat  less  cumbersome  technique  for  error  measurement 
can  be  employed  when  a  structure  is  modeled  by  elements  which 
contain  a  full  set  of  rigid  body  modes  in  their  assumed  displace¬ 
ment  fields.  Let 


be  any  rigid  body  displacement  vector  for  the  whole  structure 
(e.g.,  unit  vertical  transition  at  every  node).  Then  if  T 
is  introduced  into  the  unconstrained  force-displacement  relations 
and  the  right  hand  side  is  calculated,  there  will  result: 
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to  the  accuracy  of  the  user’s  computer.  Conversely,  solution 
of  the  constrained  force-displacement  relations: 

K„  0)  f«,7  f-KnrJ 

0  ijlrji  r„  ) 

would  result  in  #p=  /*p  to  the  accuracy  of  the  computer  if 

there  were  no  roundoff  error.  Now,  since  the  exact  can 

A  r 

be  inferred  from  T ^  merely  by  inspection,  the  above  problem 
may  be  solved  as  an  auxiliary  to  the  real  problem,  and  the  error 

measure:  , 

J  _  l*>  -Ifrl 

“  lrFi 

may  be  calculated.  The  rigid  body  mode  technique  requires  the 
saving  only  of  K  (before  constraint).  A  simple  algorithm  will 
serve  to  calculate  X“p  and  t p-  U  p,  and  the  additional  execution 
time  required  amounts  only  to  calling  each  of  BC0N,  FACTPD  (or 
FACTSD)  and  SIMULQ  once  extra.  Also,  the  lengthy  residual  force 
and  re-solution  algorithm  in  the  MAIN  program  is  avoided. 

Recent  tests  of  the  method  of  rigid  body  modes  on  a  cantilever 
beam  have  shown  that  when  dp<0..1,  dgg<<dp*  However,  when  dp>0.1 
(the  region  of  primary  interest)  dpB  performs  as  well  as  djj,. 


65 


REFERENCES 


Qrringer,  0,  "The  Effect  of  Cross  Section  Stress  Concentration 
on  the  Compressive  Strength  of  a  Unidirectional  Fiber  Composite. 
ASRL  TR  162-4  (in  preparation) 

Mack,  E.  W.,  Berg,  B.  A.,  and  Witmer,  E.  A.  "An  Improved 
Discrete-Element  Analysis  and  Program  for  the  Linear-Elastic 
Static  Analysis  of  Meridionally-Curved,  Variable-Thickness, 
Branched  Thin  Shells  of  Revolution  Subjected  to  General  Ex¬ 
ternal  Mechanical  and  Thermal  Loads,  Part  2-The  SABOR  4 
Program."  Massachusetts  Institute  of  Technology,  Aeroelastic 
and  Structures  Research  Laboratory,  ASRL  TR  146-4,  Part  2 
(also  SAMSO  TR  68-310,  Part  2),  March  1968.  (AD  840  6l4L). 

Tong,  P.  and  Plan,  T.  H.  H.  "The  Convengence  of  Finite 
Element  Method  in  Solving  Linear  Elastic.  Problems."  Int'l. 

J.  Solids  Structures,  Vol.  3,  1967.  pp.  865-879. 

Tong,  P.  "On  the  Numerical  Problems  of  the  Finite  Element 
Methods."  Study  No.  5,  Symposium  on  Computer-Aided 
Engineering,  University  of  Waterloo,  Ontario,  Canada,  1970. 

Irons,  B.  M.  and  Kan,  T.  K.  Y.  "Equation-Solving 
Algorithms  for  the  Finite-Element  Methods."  Paper  No.  V-4-2, 
University  of  Illinois  Conference  on  Numerical  Methods  in 
Structural  Analysis,  1971. 


APPENDIX  A 


C0MM0N  AREAS  DATA  REQUIREMENTS 

The  following  table  summarizes  what  information  is  expected 
by  each  FEABL  subroutine  in  the  four  control  parameter  C0MM0N 
areas  and  in  the  /DATA/  vector  C0MM0N  area. 


N.  SUBROUTINE 

X.  NAME 

C0MM0N 

AREA  X 

r 

ASEMBL 

BC0N 

FACTPD/FACTSD 

0RK 

ROTATE 

SETUP 

SIMULQ 

XTRACT 

/IO/  (Printer 

Code  KW  Only) 

/ 

/ 

/ 

/ 

/ 

/ 

/ 

■ 

/SIZE/  NET,  NDT 

/ 

/ 

/ 

/ 

/ 

/ 

■ 

/BEGIN/  Address 
index  parameters 

/ 

/ 

/ 

/ 

* 

/ 

/ 

/ 

/END/  Address 
index  parameters 

/ 

/ 

/ 

Except 

LK 

/ 

■ 

/DATA/  Vector: 

1.  Constraint  Vector 

2.  Address  Count  Vector 

3.  LNZ  Vector 

4.  Assembly  List 

5.  F/D  Vector 

6.  K  Matrix 

/ 

/ 

(c) 

/ 

■ 

/ 

/ 

/ 

/ 

/ 

II 

/ 

/ 

/ 

/ 

/ 

■ 

/ 

■ 

/ 

/ 

(a) 

/ 

m 

/ 

U 

(a) 

/ 

/ 

/ 

LDtl 

■ 

(a)  These  blocks  are  zeroed  by  0RK  before  the  first  element  is 
assembled. 

A 

(b)  Q  ~  =  Vector  of  assembled  element  equivalent  nodal  forces. 

£ 

(c)  Block  zeroed  by  SETUP 
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APPENDIX  B 


APPLYING  SHOEHORN  AND  STOPWATCH 


This  appendix  contains  data  from  which  the  user  may  estimate 
the  total  amount  of  core  storage  required  by  a  FEABL-based  pro¬ 
gram  and  the  approximate  execution  time  (CPU  time)  the  run  will 
take.  Such  estimates  will  prove  useful  aids  in  making  trade-off 
decisions.  The  data  given  in  this  appendix  is  based  on  runs  done 
on  an  IBM  370/155  using  the  F0RTPAN  G  compiler.  The  numbers  will 
vary  somewhat  from  one  machine  or  compiler  to  another. 

B.l  Estimation  of  Core  Storage  Requirement 

The  following  table  gives  the  length  of  each  FEABL  subroutine 
in  BYTES  (as  compiled  in  F0RTRAN  G)  and  words,  and  the  deck  size. 

On  IBM  360  and  370  series  hardware,  core  storage  calculations 
are  normally  done  in  terms  of  BYTES,  while  words  are  used  on  many 


Subroutine  Name 

Words 

BYTES 

No.  of  Cards 
in  Deck* 

ASEMBL 

410 

1,642 

58 

BC0N 

599 

2,396 

125 

FACTPD/FACTSD 

669 

2,674 

126 

0RK 

517 

2,070 

106 

R0TATE 

1,858 

7,430 

316 

SETUP 

726 

2,904 

94 

SIMULQ 

600 

2,398 

127 

XTRACT 

130 

522 

23 

^EABL  Software-Total 

5,509 

22,036 

975 

other  hardware  systems.  The  total  storage  requirement  for 
programs  plus  data  can  be  estimated  as  follows: 


T 


Includes  all  comment  cards 
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Item 


KWords 


KBYTES 


1.  FEABL  Software 

2.  User  programs  (a) 

3.  System  library  subprograms  (b) 

4.  /DATA/  Vector  (c) 

5.  Four  control  parameter  /C0MM0N/  areas 

5.5 

2.5 
5.1 
(L) 

22.1 

10.0 

22.0 

(4L) 

0.1 

Totals 

13.1+CL)  j  53 . 2+( 4L) 

(a)  Estimate  for  a  typical  analysis  with  a  MAIN  program 
and  two  generator  subroutines. 

(b)  Includes  library  functions  such  as  SIN,  C0S,  SQRT  and 
systems  management  subroutines. 

(c)  L=the  dimension  of  the  /DATA/  vector  in  KWords  (1000 
words) . 

B.2  Estimation  of  CPU  Time  Requirement 

The  process  of  factoring  the  master  stiffness  matrix  into 
its  triple  product: 

K  «  LDLt 


is  the  primary  time  consumer  in  any  finite  element  analysis..  Some 

study  of  the  algorithms  in  FEADL  subroutine  FACTPD/FACTSD  will 

convince  the  reader  that  the  CPU  time  consumed  is  proportional 

p 

to  NB  ,  where  N  is  the  total  number  of  unconstrained  DOF  in  the 
assembled  structure  and  B  is  the  average  semi-bandwidth  of  the 
master  stiffness  matrix.  B  may  be  calculated  approximately  from 
N  and  the  population  density  of  K  : 

n  ’  (N+1)P 


where 


P=Populatlon  Density  = 


Total  No,  of  Stored  Entries _ 

Total  Entries  in  Full  Lower  Triangle 


Experience  on  the  IBM  370/155  indicates  that  the  required  CPU 


time  is  given  approximately  by: 

time  =  NB2xlQ"^  minutes 
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APPENDIX  C 


F0RTRAN  IV  LISTING  OF  FEABL  SOFTWARE 


The  eight  subroutines  of  the  FEABL  software  system  are 
listed  in  alphabetical  order  in  this  appendix: 


Subroutine 

Page 

ASEMBL 

71 

BC0N 

73 

FACTP  D/FAC TSD 

77 

0RK 

81 

R0TATE 

84 

SETUP 

93 

SIMULQ 

96 

XTRACT 

100 

The  code  is  for  FEABL  Version  1  Release  1,  with  a  10,000-word 
/DATA/  vector.  The  following  actions  will  convert  the  program 
to  Version  2  ("On-Line")  as  explained  in  Section  2.3. 


SUBROUTINE 

Change  DIHENSI0NS 
of  REAL,  INTGR  to  2 

Delete 

Declaration 

ASEMBL 

Card  No.  0007 

■ 

Card  Nos.  0017  ana  0020 

BC0N 

0007 

0012  and  0015 

FACTP D/FACTSD 

0008 

0013  and  0016 

0RK 

0015 

0020  and  0023 

R0TATE 

0007 

0013  and  0017 

SETUP 

0007 

0013  and  0015 

SIMULQ 

0007 

0012  and  0015 

XTRACT 

0007 

0011  and  0013 
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